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V.  Introduction 
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An  investigation  of  three-dimensional  ultrasonic  mammography  is  underway.  The  goal 
of  the  research  is  improved  diagnosis  of  breast  cancer  by  quantitative,  high-resolution  three- 
dimensional  ultrasonic  imaging.  This  goal  is  being  reached  by  a  thorough  program  that 
synthesizes  recent  advances  in  tissue  modeling,  adaptive  imaging,  instrumentation,  and  sig¬ 
nal  processing.  The  final  result  of  the  research  will  be  a  major  advance  in  quantitative 
three-dimensional  ultrasonic  mammography.  Improved  resolution,  accurate  quantitative  in¬ 
formation  on  tissue  properties,  and  precise  determination  of  three-dimensional  breast  struc¬ 
ture  will  provide  crucial  new  information  for  detection,  diagnosis,  and  monitoring  of  breast 
cancer.  The  goal  of  three-dimensional  quantitative  imaging  is  currently  being  achieved  using 
novel  inverse  scattering  methods  invented  by  the  Principal  Investigator  and  coworkers.  Use  of 
full  time-domain  scattering  information  provides  images  with  high  point  resolution,  contrast 
resolution,  and  quantitative  accuracy  without  significant  artifacts.  Nonlinear  forms  of  these 
methods  provide  a  robust  approach  to  adaptive  imaging  that  is  based  on  compensation  for 
three-dimensional  scattering  from  actual  tissue  structure.  Unlike  previous  adaptive  imaging 
methods  based  on  assumptions  of  phase-screen  aberrators  and  point  scatterers,  these  meth¬ 
ods  provide  aberration  correction  ideally  suited  to  distributed  inhomogeneous  tissue  like  the 
breast.  A  unique  and  innovative  aspect  of  the  research  is  the  use  of  realistic  tissue  models 
for  ultrasonic  propagation  through  breast  tissue.  Such  models  have  been  shown  to  realisti¬ 
cally  model  wavefront  distortion  in  the  human  abdominal  wall,  but  have  not  to  date  been 
applied  to  the  human  breast.  Tissue  modeling  techniques  employ  tissue  maps  obtained  from 
specimen  cross  sections  as  well  as  from  newly  available  high-resolution  volume  photographic 
data.  Calculated  scattering  from  these  tissue  models  provides  accurate  characterization  of 
ultrasonic  propagation  within  breast  tissue  and  realistic  data  for  quantitative  imaging  al¬ 
gorithms.  The  above  studies  will  facilitate  the  application  of  breakthroughs  from  tissue 
modeling,  inverse  scattering,  and  signal  processing  to  the  critical  application  of  ultrasonic 
mammography.  Completion  of  the  proposed  research  will  make  possible  new  mammographic 
applications  of  ultrasound  that  will  provide  clinicians  with  previously  unavailable  informa¬ 
tion  and  detail.  The  end  result  will  be  a  lower-cost,  more  effective,  and  safer  modality  for 
diagnosis,  detection,  and  monitoring  of  breast  cancer. 
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VI.  Body  of  Report 

Below,  the  accomplishments  of  the  second  year  for  this  project  are  summarized  under 
the  four  categories  of  the  approved  Statement  of  Work.  Details  are  given  for  how  accom¬ 
plishments  to  date  fit  into  the  overall  research  plan.  Where  applicable,  brief  descriptions  of 
planned  research  indicate  how  the  remainder  of  the  Statement  of  Work  will  be  fulfilled. 

A.  Quantitative  Imaging  Algorithm  Development 

A  new  method  for  ultrasonic  imaging,  specifically  tailored  to  ultrasonic  mammography, 
has  also  been  developed.  This  method  allows  quantitative,  high-resolution  images  to  be 
obtained  using  direct  synthetic-aperture  processing  of  time-domain  scattered  fields.  An 
archival  manuscript  regarding  the  method  was  completed,  revised,  and  published  during  the 
second  year  of  this  project  and  is  included  as  Appendix  A  of  this  report  [1].  Additional 
results,  including  the  reconstructions  from  experimentally  measured  ultrasonic  scattering 
data,  were  presented  at  the  1999  IEEE  Ultrasonics  Symposium  in  November  1999;  a  pub¬ 
lished  manuscript  from  the  Symposium  Proceedings  [2]  is  included  as  Appendix  B.  The 
experimental  reconstructions  show  considerable  promise  for  breast  cancer  detection.  Partic¬ 
ularly  encouraging  is  the  fact  that  sub-resolvable  structures  appear  as  smoothed  variations 
rather  than  as  speckle. 

Work  has  continued  on  the  new  time-domain  imaging  method  with  implementation  non¬ 
linear  aberration  correction.  Such  correction  provides  improved  synthetic  focusing  capability 
based  on  medium  models  determined  from  quantitative  image  data.  In  one  approach,  im¬ 
proved  images  are  reconstructed  using  numerical  propagation  of  wavefields  into  estimates 
of  the  unknown  medium  [3,  4],  In  another,  much  more  efficient  approach,  the  quantitative 
sound-speed  reconstruction  obtained  by  the  time-domain  imaging  method  is  employed  to 
determine  time  shifts  that,  when  applied  to  the  measured  scattering  data,  compensate  for 
focus  aberration  caused  by  inhomogeneous  breast  tissue.  An  efficient  implementation,  in 
which  the  necessary  line  integrals  are  performed  using  a  DDA  (digital  differential  analyzer) 
method,  allows  each  iteration  to  be  performed  as  rapidly  as  the  initial  linear  reconstruction. 
Results  obtained  using  this  method  were  presented  at  the  DoD  Era  of  Hope  meeting  [5] 
(abstracts  provided  in  the  Appendices)  and  are  described  below  in  this  report. 

The  capabilities  of  the  new  time-domain  imaging  method  are  exciting  for  several  reasons. 
First,  the  images  are  both  higher  in  quality  and  more  efficiently  computed  than  conventional 
single-frequency  quantitative  images.  The  high  point  and  contrast  resolution,  as  well  as  ab¬ 
sence  of  artifacts  usually  associated  with  diffraction  tomography,  suggests  that  this  method 
will  be  very  useful  for  detection  and  characterization  of  breast  lesions.  Second,  because  of  the 
close  analogy  between  the  new  method  and  delay-and-sum  imaging,  the  new  method  could 
be  implemented  in  hardware  using  beamforming  technology  already  present  on  digital  ultra¬ 
sound  scanners.  Third,  the  quantitative  reconstructions  provided  by  the  new  method  allow 
aberration  correction  to  be  implemented  much  more  robustly  than  possible  in  conventional 
pulse-echo  imaging.  The  new  time-domain  method  also  has  the  capability  to  incorporate 
other  imaging  techniques  ( e.g .,  time-gain  compensation  and  harmonic  imaging)  currently 
used  in  clinical  and  experimental  B-scan  systems. 
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B.  Tissue  Modeling 

Progress  toward  improved  scattering  models  for  ultrasound-breast  tissue  interaction  has 
been  made  in  several  studies. 

Work  on  improved  scattering  models  for  ultrasound-tissue  interaction  continued.  A  paper 
on  simulation  of  ultrasonic  propagation  through  cross-sectional  models  of  chest  wall  tissue 
was  revised  and  published  in  JASA  [6]  and  is  included  here  as  Appendix  C.  This  work  expands 
upon  previous  models  by  including  tissue-dependent  absorption  effects  and  by  analysis  of 
the  frequency  dependence  of  ultrasonic  wavefront  distortion. 

Work  on  fc-space  methods  for  measurements  of  ultrasonic  scattering  also  continued,  and 
significant  new  progress  was  made.  A  collaborative  project  with  Weidlinger  Associates  and 
the  University  of  Rochester  provided  detailed  quantitative  comparisons  between  a  new  k- 
space  method  and  a  state-of-the-art  pseudospectral  solver.  This  study  showed  that  the 
k- space  method  has  strong  advantages  for  large-scale  simulations  of  propagation  through 
soft  tissue,  which  is  precisely  the  simulation  problem  of  interest  for  the  USAMRMC-funded 
breast  cancer  study.  The  comparison  was  presented  at  the  1999  IEEE  Ultrasonics  Symposium 
and  has  been  published  in  the  proceedings  of  that  symposium;  a  copy  of  that  publication  is 
included  here  as  Appendix  D  [7]. 

Knowledge  gained  from  the  abovementioned  comparison  project  led  to  a  number  of  im¬ 
provements  in  the  A;-space  method,  including  improved  computational  efficiency,  more  ac¬ 
curate  interpolation  of  simulated  pressure  signals,  and  more  effective  methods  to  compute 
scattering  from  media  including  discontinuities.  These  improvements  were  reported  in  a  pre¬ 
sentation  to  the  Acoustical  Society  of  America  [4],  (abstract  provided  in  the  Appendices). 
These  comparisons,  as  well  as  new  analysis  that  explains  the  remarkable  stability  and  accu¬ 
racy  of  the  /e-space  method,  greatly  improved  the  extensively  revised  version  of  a  manuscript 
accepted  for  publication  [8]  that  is  provided  here  as  Appendix  E. 

Now  that  a  robust,  efficient  three-dimensional  method  for  computation  of  ultrasonic 
propagation  is  available,  work  on  mapping  of  breast  tissue  has  begun  in  earnest.  In  collabo¬ 
ration  with  colleagues  at  the  University  of  Rochester,  cross-sectional  breast  tissue  specimens 
have  been  sectioned  and  stained  for  high-resolution  segmentation  by  tissue  type.  One  such 
tissue  cross  section  is  shown,  together  with  a  computed  ultrasonic  pulse  propagating  within 
the  tissue,  in  Fig.  1.  The  ultrasonic  propagation  was  computed  using  the  ft-space  method, 
which  provides  much  greater  efficiency  and  accuracy  than  the  finite-difference  method  used 
in  previous  simulation  studies  of  propagation  through  tissue  [6,  9,  10]. 

During  the  second  year  of  the  project,  photographic  image  data  from  the  Visible  Woman 
data  set  [11]  was  also  analyzed  to  obtain  quantitative  maps  of  breast  tissue  for  simulation 
studies  and  testing  of  imaging  algorithms.  A  method  in  which  hue,  saturation,  and  value  are 
mapped  to  tissue  type  was  applied;  this  method  also  incorporates  nonlinear  processing  to 
enforce  uniformity  among  multiple  layers  in  three-dimensional  breast  models.  An  example 
two-dimensional  simulation,  using  the  A;-space  method  and  a  Visible  Woman  breast  model, 
is  shown  in  Fig.  2. 

During  the  third  year  of  the  project,  the  new  automatic  segmentation  method  will  be 
applied  to  three-dimensional  photographic  data  from  the  Visible  Woman  data  set  to  ob¬ 
tain  definitive  three-dimensional  models  of  breast  tissue  structure  for  studies  of  ultrasonic 
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Figure  1:  Ultrasonic  propagation  of  a  2.5  MHz  pulse  through  a  two-dimensional  breast  cross- 
sectional  model  obtained  from  a  segmented,  stained  breast  tissue  cross  section.  An  area  of 
167  x  50  mm2  is  shown. 


Figure  2:  Three  frames  showing  ultrasonic  propagation  of  a  1  MHz  pulse  through  a  two- 
dimensional  breast  cross-sectional  model  obtained  from  the  Visible  Woman  data  set.  Each 
panel  shows  an  area  of  97  x  69  mm2. 


propagation  and  imaging. 

C.  Quantitative  Imaging  Algorithm  Implementation 

Quantitative  reconstructions  performed  using  the  time-domain  diffraction  tomography 
method  are  shown  in  Refs.  [1]  (Appendix  A)  and  [2]  (Appendix  B)  for  simulated  two- 
dimensional  and  three-dimensional  scattering  data.  These  results  are  promising  for  ultrasonic 
mammography.  As  discussed  below  in  section  D  (Evaluation  and  Comparison  of  Results), 
time-domain  quantitative  images  show  parametric  accuracy,  high  resolution,  and  few  arti¬ 
facts.  Effects  of  limited  scattering  data,  shown  in  Fig.  3  of  Appendix  A  for  the  2D  case, 
indicate  that  accurate  images  can  be  obtained  without  the  necessity  of  apertures  that  entirely 
enclose  the  breast. 

Reconstructions  employing  experimentally  measured  scattering  data  have  also  been  per¬ 
formed  in  collaboration  with  colleagues  at  the  University  of  Rochester.  Tissue-mimicking 
phantoms  composed  of  agar  gel  with  glass  spheres  have  been  constructed  and  their  scat¬ 
tering  has  been  measured  using  a  2048-element  ring  transducer.  Reconstructions  using  the 
new  time-domain  diffraction  tomography  method  are  shown  in  Ref.  [2]  (Appendix  B).  These 
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Figure  3:  Aberration  correction  applied  to  time-domain  quantitative  imaging  of  a  breast 
tissue  model.  From  left  to  right:  unaberrated  image;  first  iteration;  second  iteration.  Each 
panel  shows  the  reconstructed  sound  speed  as  a  linear  grayscale  image,  with  black  indicating 
the  maximum  sound  speed  ( e.g .,  connective  tissue)  and  white  indicating  minimum  speed 
{e.g.,  fat). 

reconstructions  show  great  promise  for  practical  implementation  of  the  new  time-domain 
diffraction  tomography  method  for  breast  imaging  in  vivo. 

An  anthropomorphic  breast-mimicking  phantom,  developed  in  collaboration  between 
General  Electric,  University  of  Wisconsin,  and  University  of  Rochester,  is  also  available 
for  testing.  Time-domain  scattering  data,  to  be  taken  with  the  2048-element  University 
of  Rochester  ring  transducer  [12]  will  be  used  as  input  for  testing  of  the  new  time-domain 
quantitative  imaging  method.  Both  2D  and  3D  reconstructions  will  be  performed. 

As  described  above,  a  nonlinear,  aberration-corrected  version  of  the  new  time-domain 
inverse  scattering  method  has  been  implemented  and  tested  with  simulated  data  (from  the 
breast  tissue  models  described  above  in  section  B,  Tissue  Modeling).  Preliminary  results, 
as  seen  in  Figure  3,  indicate  that  the  aberration  correction  scheme  will  allow  accurate  ul¬ 
trasonic  images  of  whole  breasts  to  be  obtained.  During  the  third  year  of  the  project, 
aberration-corrected  time-domain  diffraction  tomography  will  also  be  applied  to  imaging  of 
large  phantoms,  such  as  the  GE  breast  phantom,  from  experimentally  measured  scattering 
data. 

D.  Evaluation  and  Comparison  of  Results 

Analysis  of  the  time-domain  imaging  results  for  simulated  data  is  given  in  Ref.  [1]  (Ap¬ 
pendix  A).  These  results  indicate  that  high  quantitative  accuracy  can  be  achieved.  Compu¬ 
tations  of  point-spread  functions  (Fig.  2  in  Appendix  A)  show  that  the  time-domain  method 
yields  higher  point  and  contrast  resolution  than  single-frequency  diffraction  tomography.  For 
the  3D  case,  the  level  of  the  first  sidelobe  is  reduced  by  13  dB,  while  the  second  sidelobe  is 
reduced  by  18  dB.  Because  of  the  broadband  scattering  information  employed,  the  width  of 
the  main  lobe  indicates  point  resolution  of  features  smaller  than  one-half  wavelength  at  the 
center  frequency. 

Evaluation  of  quantitative  accuracy  (Fig.  4  in  Appendix  A)  shows  that  the  time-domain 
diffraction  tomography  method  provides  parametric  accuracy  similar  to  established  single- 
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frequency  methods  while  showing  much  more  immunity  from  artifacts.  The  new  aberration 
correction  method,  described  above  and  illustrated  in  Fig.  1,  provides  high  quantitative 
accuracy  even  for  imaging  of  large  inhomogeneities  for  which  the  Born  approximation  does 
not  hold. 

Two-dimensional  and  three-dimensional  reconstructions  performed  using  the  new  time- 
domain  quantitative  imaging  method  will  be  directly  compared  to  analogous  images  from 
other  modalities  such  as  x-ray  computed  tomography  and  conventional  B-scan  ultrasound. 
Consultants  Dr.  Jonathan  Meilstrup  and  Dr.  Claudia  Kasales  will  help  with  this  portion  of 
evaluation  and  comparison  of  results.  This  evaluation  task  will  take  place  during  the  third 
year  of  the  project. 

Quantitative  evaluation  of  the  k- space  method  for  simulation  of  ultrasonic  propagation 
has  been  performed  and  is  reported  in  Refs.  [7]  (Appendix  D)  and  [8]  (Appendix  E).  These 
results  show  that  the  /c-space  method  provides  much  higher  accuracy  than  finite-difference 
or  pseudospectral  methods.  A  particularly  striking  result  is  that  simulated  ultrasonic  waves 
can  propagate  many  thousands  of  wavelengths  without  distortion  [4]. 

Regarding  the  task  of  publishing  results  in  archival  journals,  the  second  year  of  the  project 
has  resulted  in  two  published  articles,  one  accepted  article,  and  two  articles  in  a  peer-reviewed 
conference  proceedings  volume.  Further  publications  planned  to  be  prepared  during  the  next 
year  will  be  based  on  abberation-corrected  time-domain  diffraction  tomography,  modeling  of 
breast  tissue  based  on  the  Visible  Woman  data  set,  and  scattering  computations  employing 
realistic  breast  tissue  models. 
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VII.  Key  Research  Accomplishments 

The  key  research  accomplishments  to  date  in  this  project  can  be  summarized  as  follows: 

•  Development  of  a  new  time-domain  quantitative  imaging  method  designed  specifically 
for  ultrasonic  mammography. 

•  Numerical  implementation  and  testing  of  the  new  time-domain  imaging  method,  show¬ 
ing  that  this  method  provides  high  accuracy  with  greater  efficiency  than  previous  in¬ 
verse  scattering  methods. 

•  Testing  of  the  time-domain  imaging  method  using  scattering  data  measured  by  a  2048- 
element  ring  transducer. 

•  Implementation  of  abberation  correction  to  the  time-domain  quantitative  imaging 
method. 

•  Implementation  of  a  tissue-dependent  absorption  model  for  simulations  of  scattering 
and  propagation. 

•  Characterization  of  the  frequency  dependence  of  ultrasonic  scattering  from  human  soft 
tissues. 

•  Exact  computation  of  time-domain  scattering  from  simple  objects  for  testing  of  quan¬ 
titative  imaging  methods. 

•  Implementation  and  testing  of  a  new  fc-space  method  for  computation  of  scattering, 
indicating  that  the  method  is  accurate  and  extremely  efficient,  and  therefore  ideal  for 
the  proposed  3D  computations  of  scattering  from  breast  tissue. 

•  Extension  of  the  new  /c-space  method  to  include  tissue-dependent  absorption,  absorbing 
boundary  layers,  and  three-dimensional  scattering. 

•  Development  of  realistic  breast  tissue  models  using  stained  tissue  cross  sections  and 
photographic  data  from  the  Visible  Woman  data  set. 

•  Use  of  the  k- space  method  to  compute  time-domain  scattering  data  from  realistic  breast 
tissue  models,  for  analysis  of  wavefront  distortion  and  testing  of  quantitative  imaging 
algorithms. 

•  Quantitative  analysis  of  inverse  scattering  results,  indicating  that  time-domain  recon¬ 
structions  provide  much  higher  point  resolution,  contrast  resolution,  and  freedom  from 
artifacts  than  single-frequency  reconstructions. 
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VIII.  Reportable  Outcomes 

Reportable  outcomes  for  the  second  year  of  this  project  research  have  included  two  pa¬ 
pers  published  in  archival  journals  [1,  6],  one  additional  accepted  paper  [8],  two  papers 
published  in  the  2000  IEEE  Ultrasonics  Symposium  Proceedings  [2,  7],  and  two  published 
abstracts  presented  at  scientific  meetings  [5,  4].  All  of  these  publications,  as  well  as  a  current 
Curriculum  Vitae  for  the  Principal  Investigator,  are  included  below  in  the  Appendices. 

Accomplishments  performed  under  the  USAMRMC-funded  project  have  facilitated  ap¬ 
plications  for  other  funding  support.  The  proposal  “High-Resolution  Breast  Tissue  Mapping 
using  Pulse-Echo  Ultrasonography”  (T.  Douglas  Mast,  Principal  Investigator)  was  submit¬ 
ted  to  the  Concept  Award  program  of  the  USAMRMC.  This  study  was  not  funded.  An 
additional  proposal,  “Optimized  Intracavitary  Ultrasound  Array  for  Uniform  Hyperther¬ 
mia  Treatment  of  Prostate  Cancer”  (Nadine  B.  Smith,  Principal  Investigator;  T.  Douglas 
Mast,  Co-Principal  Investigator;  funding  notification  pending)  was  submitted  to  the  DoD 
Prostate  Cancer  Research  Program.  Although  not  a  direct  extension  of  the  present  ultrasonic 
mammography  project,  the  proposed  project  would  leverage  the  advanced  tissue  modeling 
techniques  developed  for  the  current  USAMRMC-funded  research. 

The  Principal  Investigator  of  this  project  has  been  appointed  to  the  position  of  Assistant 
Professor  in  the  Pennsylvania  State  University  Graduate  Program  in  Acoustics.  During  the 
Spring  of  2000,  he  taught  a  well-received  upper-level  graduate  course  on  acoustic  scattering. 
This  course  included  a  comprehensive  treatment  of  ultrasonic  scattering  by  human  tissues, 
and  was  greatly  enriched  by  the  Principal  Investigator’s  research  performed  for  the  present 
USAMRMC-funded  project. 

A  student,  James  F.  Kelly,  has  recently  been  added  to  the  project  team.  Mr.  Kelly  is 
not  directly  funded  by  the  USAMRMC,  but  instead  is  supported  by  the  Mathematics  Honor 
Student  program  of  the  Applied  Research  Laboratory.  He  is  currently  exploring  aspects  of  the 
time-domain  inverse  scattering  problem,  including  deconvolution  of  scattered  wavefields,  that 
extend  the  present  Statement  of  Work  (some  of  this  includes  research  into  ideas  discussed 
in  the  abovementioned  Concept  Award  proposal).  If  successful,  these  exploratory  studies 
will  allow  the  imaging  methods  developed  under  the  USAMRMC-funded  research  to  achieve 
practical  application  even  more  rapidly. 
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IX.  Conclusions 


The  second  year  of  this  USAMRMC-funded  project  has  yielded  considerable  progress  to¬ 
ward  the  goal  of  improved  diagnosis  of  breast  cancer  by  three-dimensional  ultrasonic  imaging. 
Several  breakthroughs  have  been  made  which  will  provide  a  solid  foundation  for  the  con¬ 
tinued  tissue  modeling  and  quantitative  imaging  research  planned  for  the  last  year  of  the 
project. 

Breakthroughs  in  the  area  of  quantitative  imaging  have  included  successful  imaging  of 
tissue-mimicking  phantoms  using  measured  scattering  data  as  well  as  implementation  of  a 
new  aberration  correction  method  for  quantitative  ultrasonic  mammography.  The  present 
method  is  potentially  of  very  great  importance  for  breast  cancer  diagnosis  for  several  rea¬ 
sons:  (1)  images  have  higher  quality  than  that  achievable  by  conventional  inverse  scattering 
methods  or  by  current  ultrasound  scanners,  (2)  tissue  parameters  are  computed  and  quanti¬ 
tatively  imaged  with  high  accuracy,  and  (3)  the  close  analogy  between  the  new  method  and 
conventional  synthetic-aperture  imaging  will  allow  rapid  implementation  of  the  new  method 
on  hardware  similar  to  currently  used  beamformers.  The  additional  improvement  of  aberra¬ 
tion  correction  increases  the  value  of  this  method  even  more,  because  the  strong  scattering 
inherent  to  breast  tissue  [13]  is  an  important  limiting  factor  to  existing  ultrasonic  imaging 
methods. 

In  the  area  of  breast  tissue  modeling,  simulation  methods  have  been  developed  that  will 
allow  efficient  computation  of  scattering  from  fully  three-dimensional  models  of  breast  tissue. 
A  new  A;-space  method  has  been  implemented,  rigorously  tested,  and  extended  to  include  ab¬ 
sorbing  boundary  conditions  as  well  as  tissue-dependent  ultrasonic  absorption.  The  method 
has  been  shown  to  provide  the  accuracy  and  efficiency  needed  for  large-scale  computations 
of  ultrasonic  propagation  within  the  human  breast.  A  three-dimensional  implementation  of 
this  method  is  already  in  place.  Breast  tissue  models  have  been  developed  using  both  stained 
cross  sections  of  human  breast  tissue  and  photographic  image  data  from  the  Visible  Woman 
project  of  the  National  Library  of  Medicine. 

Future  work  based  on  these  breakthroughs  will  follow  the  approved  Statement  of  Work. 
Plans  for  the  final  year  of  research  include  detailed  simulation  of  ultrasound  interaction  with 
breast  tissue  for  two-  and  three-dimensional  tissue  models.  Definitive  aberration-corrected 
reconstructions  will  be  computed  both  for  measured  scattering  data  from  breast-mimicking 
phantoms  and  for  simulated  scattering  data  using  detailed  breast  models.  These  results 
will  be  carefully  evaluated  for  accuracy  and  clinical  utility,  in  collaboration  with  the  named 
clinical  consultants. 

The  final  outcome  of  the  successfully  completed  project  will  be  a  novel  method  for  early 
detection,  characterization,  and  treatment  monitoring  of  breast  cancer  lesions.  Results  to 
date  indicate  that  the  final  method  will  provide  image  quality  and  diagnostic  information 
greatly  superior  to  current  2D  and  3D  ultrasonic  mammography  methods  The  finally  re¬ 
sulting  method  is  expected  to  be  competitive  with  magnetic  resonance  imaging  and  x-ray 
computed  tomography  as  a  tool  for  breast  cancer  diagnosis,  while  maintaining  inherent  ad¬ 
vantages  of  ultrasound  such  as  lower  cost,  ability  to  characterize  cystic  and  solid  lesions,  and 
safe  nonionizing  radiation. 
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A  quantitative  ultrasonic  imaging  method  employing  time-domain  scattering  data  is  presented.  This 
method  provides  tomographic  images  of  medium  properties  such  as  the  sound  speed  contrast;  these 
images  are  equivalent  to  multiple-frequency  filtered-backpropagation  reconstructions  using  all 
frequencies  within  the  bandwidth  of  the  incident  pulse  employed.  However,  image  synthesis  is 
performed  directly  in  the  time  domain  using  coherent  combination  of  far-field  scattered  pressure 
waveforms,  delayed  and  summed  to  numerically  focus  on  the  unknown  medium.  The  time-domain 
method  is  more  efficient  than  multiple-frequency  diffraction  tomography  methods,  and  can,  in  some 
cases,  be  more  efficient  than  single-frequency  diffraction  tomography.  Example  reconstructions, 
obtained  using  synthetic  data  for  two-  and  three-dimensional  scattering  of  wideband  pulses,  show 
that  the  time-domain  reconstruction  method  provides  image  quality  superior  to  single-frequency 
reconstructions  for  objects  of  size  and  contrast  relevant  to  medical  imaging  problems  such  as 
ultrasonic  mammography.  The  present  method  is  closely  related  to  existing  synthetic-aperture 
imaging  methods  such  as  those  employed  in  clinical  ultrasound  scanners.  Thus,  the  new  method  can 
be  extended  to  incorporate  available  image-enhancement  techniques  such  as  time-gain 
compensation  to  correct  for  medium  absorption  and  aberration  correction  methods  to  reduce  error 
associated  with  weak  scattering  approximations.  ©  1999  Acoustical  Society  of  America. 
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PACS  numbers:  43.20.Fn,  43.60.Rw,  43.80. Vj,  43.20.Px  [ANN] 


INTRODUCTION 

Quantitative  imaging  of  tissue  properties  is  a  potentially 
useful  technique  for  diagnosis  of  cancer  and  other  pathologi¬ 
cal  conditions.  Inverse  scattering  methods  such  as  diffraction 
tomography  can  provide  quantitative  reconstruction  of  tissue 
properties  including  sound  speed,  density,  and  absorption. 
However,  although  previous  inverse  scattering  methods  have 
achieved  high  resolution  and  quantitative  accuracy,  such 
methods  have  not  yet  been  incorporated  into  commercially 
successful  medical  ultrasound  imaging  systems. 

Current  inverse  scattering  methods  are  lacking  in  several 
respects  with  respect  to  conventional  B-scan  and  synthetic 
aperture  imaging  techniques.  Previous  methods  of  diffraction 
tomography,  including  methods  based  on  the  Bom  and  Ry- 
tov  approximations,1,2  and  higher-order  nonlinear 
approaches,3,4  have  usually  been  based  on  single-frequency 
scattering,  while  current  diagnostic  ultrasound  scanners  em¬ 
ploy  wideband  time-domain  signals.  The  use  of  wideband 
information  in  image  reconstruction  is  known  to  provide  in¬ 
creased  point  and  contrast  resolution,5,6  both  of  which  are 

A  c  7  o 

important  for  medical  diagnosis.  ’  ’ 

Several  approaches  have  been  used  to  incorporate  wide¬ 
band  scattering  information  into  quantitative  ultrasonic  im¬ 
aging.  One  group  of  methods  employs  time-domain  tomog¬ 
raphy  based  on  Radon-transform  relationships  that  hold 
(under  the  assumption  of  weak  scattering)  between  scattered 
acoustic  fields  and  the  reflectivity  or  scattering  strength  of 
the  medium.  Pioneering  work  in  this  area9'10  employed  mea¬ 
surements  of  reflectivity  in  pulse-echo  mode,  while  later 
studies  have  incorporated  aberration  correction11'12  and 
multiple-angle  scattering  measurements.13,14  A  limitation  of 


these  methods,  however,  is  that  the  Radon  transform  rela¬ 
tionship  strictly  holds  only  when  the  medium  is  insonified  by 
an  impulsive  (infinite  bandwidth)  wave.  When  pulses  of  fi¬ 
nite  bandwidth  are  employed,  image  quality  can  degrade 
significantly.15 

A  number  of  linear  and  nonlinear  diffraction  tomogra¬ 
phy  methods  have  been  implemented  using  scattering  data 
for  a  number  of  discrete  frequencies  (e.g.,  Refs.  16-19).  Al¬ 
though  use  of  multiple-frequency  data  provides  improve¬ 
ments  in  image  quality,  computational  requirements  for 
multiple-frequency  imaging  are  typically  large  because  the 
computational  cost  is  proportional  to  the  number  of  frequen¬ 
cies  employed.  To  achieve  image  quality  competitive  with 
present  diagnostic  scanners,  together  with  quantitative  imag¬ 
ing  of  tissue  properties,  present  frequency-domain  methods 
may  require  solution  of  the  inverse  scattering  problem  for 
many  frequencies  within  the  bandwidth  of  the  transducer 
employed.  This  approach  thus  demands  a  high  computational 
cost,  so  that  high-quality  real-time  imaging  may  not  be  pres¬ 
ently  feasible  using  current  frequency-domain  inverse  scat¬ 
tering  methods. 

Very  few  previous  workers  have  investigated  direct  use 
of  time-domain  waveform  data  for  inverse  scattering  meth¬ 
ods  analogous  to  frequency-domain  diffraction  tomography. 
Several  methods20  21  have  used  frequency  decomposition  of 
scattered  pulses  to  construct  a  wideband  estimate  of  the  spa¬ 
tial  Fourier  transform  of  an  unknown  medium;  after  appro¬ 
priate  averaging  and  interpolation,  this  transform  can  be  in¬ 
verted  to  obtain  a  wideband  Born  reconstruction  of  the 
medium.  A  study  reported  in  Ref.  22  has  showed  that  broad¬ 
band  synthetic  aperture  imaging  using  linear  arrays  is  closely 


3061  J.  Acoust.  Soc.  Am.  106  (6),  December  1999  0001  -4966/99/1 06(6)73061/1 1/$ 15.00  ©  1999  Acoustical  Society  of  America  3061 


related  to  inverse  scattering  using  filtered  backpropagation. 
A  related  method,  suggested  in  Ref.  23,  provides  a  time- 
domain  reconstruction  algorithm  that  employs  filtered  back- 
propagation  of  scattered  waveforms  measured  on  a  circular 
boundary.  However,  the  time  domain  reconstruction  formula 
of  Ref.  23  yields  reconstructions  that  are  less  general  than 
multiple-frequency  reconstructions  obtained  using  the  same 
signal  bandwidth. 

Another  approach,  related  both  to  multiple-frequency 
methods  and  direct  time-domain  methods,  has  recently  been 
presented.24  This  work  extends  the  eigenfunction  method  of 
Ref.  19  to  use  the  full  bandwidth  of  the  incident  pulse  wave¬ 
form.  In  the  extended  method,  eigenfunctions  and  eigenval¬ 
ues  of  a  scattering  operator  are  computed  to  obtain  a 
frequency-dependent  representation  of  the  scattering  me¬ 
dium.  Fourier  synthesis  is  then  applied  to  obtain  a  time- 
dependent  estimate  of  the  medium.  A  cross-correlation  op¬ 
eration  removes  the  time  dependence  of  the  estimate  as  well 
as  its  dependence  on  the  waveform  employed. 

The  present  paper  offers  a  new  approach  to  wideband 
quantitative  imaging:  a  time-domain  inverse  scattering 
method  that  overcomes  some  of  the  limitations  of  previous 
frequency-domain  and  time-domain  quantitative  imaging 
methods.  The  new  method  provides  tomographic  reconstruc¬ 
tions  of  unknown  scattering  media  using  the  entire  available 
bandwidth  of  the  signals  employed.  Reconstructions  are  per¬ 
formed  using  scattering  data  measured  on  a  surface  sur¬ 
rounding  the  region  of  interest,  so  that  the  method  is  well 
suited  to  ultrasonic  mammography.  The  reconstruction  algo¬ 
rithm  is  derived  as  a  simple  delay-and-sum  formula  similar 
to  synthetic-aperture  algorithms  employed  in  conventional 
clinical  scanners.  However,  unlike  current  clinical  scanners, 
the  present  method  can  provide  quantitative  images  of  tissue 
properties  such  as  the  spatially  dependent  sound  speed.  Re¬ 
constructions  obtained  in  this  manner  are  equivalent  to  re¬ 
constructions  obtained  by  combining  conventional 
frequency-domain  diffraction  tomography  reconstructions 
for  all  frequencies  within  the  signal  bandwidth  of  interest. 
The  current  method,  however,  can  be  even  more  efficient 
than  single-frequency  diffraction  tomography.  The  method  is 
applicable  both  to  two-dimensional  and  three-dimensional 
image  reconstruction.  The  direct  time-domain  nature  of  the 
reconstruction  algorithm  allows  straightforward  incorpora¬ 
tion  of  depth-  and  frequency-dependent  amplitude  correction 
to  compensate  for  medium  absorption  as  well  as  aberration 
correction  methods  to  overcome  limitations  of  the  Born  ap¬ 
proximation. 

I.  THEORY 

A.  The  time-domain  reconstruction  algorithm 

An  inverse  scattering  algorithm,  applicable  to  quantita¬ 
tive  imaging  of  tissue  and  other  inhomogeneous  media,  is 
derived  below.  For  simplicity  of  derivation,  the  medium  is 
modeled  as  a  fluid  medium  defined  by  the  sound  speed  con¬ 
trast  function 


y(r)  = 


c(r)2 


-1, 


(1) 


FIG.  1.  Scattering  configuration.  An  incident  pressure  pulse  f(t—  a •  rlc )  is 
scattered  by  an  inhomogeneous  medium  and  the  time-domain  scattered  pres¬ 
sure  ps(0,a^t)  is  measured  at  a  radius  R  in  the  far  field. 


where  c0  is  a  background  sound  speed  and  c(r)  is  the  spa¬ 
tially  dependent  sound  speed  defined  at  all  points  r.  For  the 
scope  of  the  initial  derivation,  the  medium  is  assumed  to 
have  constant  density,  no  absorption,  and  weak  scattering 
characteristics;  extensions  to  the  reconstruction  algorithm 
that  overcome  these  limiting  assumptions  are  discussed  in 
the  following  section. 

For  the  model  of  the  scattering  medium  represented  by 
Eq.  (1),  the  time-domain  scattered  acoustic  pressure  ps(r,t) 
obeys  the  wave  equation25 

1  d2ps{r,t)  _  y(r)  d2p(r,t) 

9  9  9  9  >  \"/ 

Cq  Sr  Cq  Sr 

where  p(r,t )  is  the  total  acoustic  pressure  in  the  medium. 

The  scattering  configuration  considered  here  is  sketched 
in  Fig.  1.  The  medium  is  subjected  to  a  pulsatile  plane  wave 
propagating  in  the  direction  of  the  unit  vector  a , 

pinc(r,a,t)=f(t-r-a/c0),  (3) 

where  /  is  the  time-domain  waveform  and  c0  is  the  back¬ 
ground  sound  speed.  The  scattered  wavefield  ps(0,a,t)  is 
measured  at  a  fixed  radius  R  in  the  far  field,  where  0  corre¬ 
sponds  to  the  direction  unit  vector  of  a  receiving  transducer 
element.  (Alternatively,  if  scattering  measurements  are  made 
in  the  near  field,  the  far-field  acoustic  pressure  can  be  com¬ 
puted  using  exact  transforms  that  represent  propagation 
through  a  homogeneous  medium.16) 

A  general  time-domain  solution  for  the  wave  equation 
(2),  valid  for  two-dimensional  (2D)  or  three-dimensional 
(3D)  scattering,  is  then 


V2/»,(r,0- 


ps(0,a,t)=  j  ps(0,a,a))e  loJtdco, 

J  -oo 


(4) 


where  ps(0,a,a))  is  a  single  frequency  component  of  the 
scattered  wavefield, 

ps(0,a,(o)  =  J  ^  ps(0,a,t)eia”dt,  (5) 

given  exactly  by25 

ps(0,a,(o)  =  k2f(w)  J  Go(R0- r0,o>) 


Xy(r0)p(r,a,a))dV0.  (6) 

In  Eq.  (6),  k  is  the  wave  number  cd/c0  and  p(r0,a,a))  is  the 
total  acoustic  pressure  associated  with  the  unit-amplitude  in¬ 
cident  plane  wave  elka'r°.  The  integral  in  Eq.  (6)  is  taken 
over  the  entire  support  of  y  in  R2  for  2D  scattering  or  in  R3 
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(10) 


foV  3D  scattering.  The  free-space  Green’s  function,  repre¬ 
sented  by  G0  in  Eq.  (6),  is26 

G0(r,o>)  =  -H^fcr)  for  2D  scattering 


and  (7) 

eikr 

G0(  r,  a> )  =  for  3D  scattering, 

where  Hq!)  is  the  zeroth-order  Hankel  function  of  the  first 
kind  and  r  is  the  magnitude  of  the  vector  r. 

The  far-field  scattered  pressure,  when  specified  for  all 
incident-wave  directions  a,  measurement  directions  0 \  and 
times  t,  comprises  the  data  set  to  be  used  for  reconstruction 
of  the  unknown  medium.  The  inverse  scattering  problem, 
specified  by  Eq.  (6)  for  a  single  frequency  component,  is  to 
reconstruct  the  unknown  medium  contrast  y(  r)  using  the 
measured  data  ps(0,a,(o). 

The  starting  point  for  the  present  time-domain  inverse 
scattering  method  is  conventional  single-frequency  diffrac¬ 
tion  tomography.  Under  the  assumption  of  weak  scattering, 
one  can  make  the  Bom  approximation,  in  which  the  total 
pressure  p(a,(o )  in  Eq.  (6)  is  replaced  by  the  plane  wave 
eikr  a.  For  scattering  measurements  made  at  a  radius  R  in  the 
far  field,  the  linearized  inverse  problem  of  Eq.  (6)  can  be 
then  solved  for  any  frequency  component  using  filtered 
backpropagation, 2,16,27  i.e., 

p(a>)e~ikR  f  f 

yB(r,<o)  = — r— —  $>(0,a) 

f((o)  J  J 


where 

£(«»)=- 


X ps(  0 ,  a ,  <a)eikie~a)  TdSadSe 


<£(0,a)  =  |sin(0—  a)|  in  2D, 

and 


(8) 


(9) 


/2(o))= — <t>(0,a)  =  \0- a\  in  3D. 

4  7i 

Each  surface  integral  in  Eq.  (8)  is  performed  over  the  entire 
measurement  circle  for  the  2D  case  and  over  the  entire  mea¬ 
surement  sphere  for  the  3D  case.  Equation  (8)  provides  an 
exact  solution  to  the  linearized  inverse  scattering  problem  for 
a  single  frequency  component  of  the  scattered  wavefield 
ps(0,a,t).  The  resulting  reconstruction,  yB(r,(o),  has  spatial 
frequency  content  limited  by  the  “Ewald  sphere”  of  radius 
2  k  in  wavespace.1 

To  improve  upon  the  single-frequency  formulas  speci¬ 
fied  by  Eq.  (8),  one  can  extend  the  spatial-frequency  content 
of  reconstructions  by  exploiting  wideband  scattering  infor¬ 
mation.  The  method  outlined  here  synthesizes  a  “multiple- 
frequency”  reconstruction  yM{ r)  by  formally  integrating 
single- frequency  reconstructions  y5(r,o>)  over  a  range  of 
frequencies  a) .  A  generalized  formula  for  this  approach  can 
be  written 


7u(  r)  = 


fog(v)da) 


where  g((o )  is  an  appropriate  frequency-dependent  weight¬ 
ing  function.  In  practice,  the  weighting  function  g(a>)  is  cho¬ 
sen  to  be  bandlimited  because  (for  a  given  set  of  physical 
scattering  measurements)  the  frequency-dependent  contrast 
yB( r,o>)  can  only  be  reliably  reconstructed  for  a  finite  range 
of  frequencies  o)  associated  with  the  spectra  of  the  incident 
waves  employed.  Thus,  the  integrands  in  Eq.  (10)  are  non¬ 
zero  only  over  the  support  of  g(w)  and  the  corresponding 
integrals  are  finite. 

Using  Eq.  (8),  and  making  the  definition 


(ID 


Eq.  (10)  can  be  written  in  the  form 


r*f(r)  = 


2 

N 


fi.(a>)e-‘kR 

/(*>) 


<*>(*«) 


X  ps(O,a,w)eik<0~a)rdSadS  gd(a.  (12) 

If  the  frequency  weight  g((o)  is  now  specified  to  incor¬ 
porate  the  incident-pulse  spectrum /(w)  and  to  compensate 
for  the  frequency-  and  dimension-dependent  coefficient 


£(&>)= 


/(<*>) 

£(«)’ 


Eq.  (12)  reduces  to  the  form 


(13) 


rA/(r)=^J  J  <p(0,a)f°  ps(0,a,a> ) 


Xc-/*[*+(«-0)  r ld(odSadS0.  (14) 

The  choice  of  frequency  weight  from  Eq.  (13)  allows  the 
multiple-frequency  reconstruction  formula  of  Eq.  (12)  to  be 
greatly  simplified.  Specifically,  the  inner  integral  of  Eq.  (14) 
resembles  a  weighted  inverse  Fourier  transform  of  the 
frequency-domain  scattered  field  p(0,a,a>).  To  obtain  an 
explicit  time-domain  expression  for  7M(r)>  Eq.  (14)  can  be 
rewritten  using  the  definition  of  ps(0,a,(o)  from  Eq.  (5)  to 
yield 


0,a,Rlco+ 


(«-<)•  r 

co  i 


dSadS0 ,  (15) 


where  L  denotes  the  linear  operator 


LOW] 


-2f 


if/((x))e  ltotda) 


(16) 


and  $/(<*))  is  the  Fourier  transform  of  ijj{t)  using  the  defini¬ 
tion  from  Eq.  (5). 

Using  the  conjugate  symmetry  of  <A(co)  [i.e., 
lj/(  0 ;  a, ay)  =  {//*(  0 ;  a,  —  o) )  for  any  real  i/s(t)] ,  the  real  part  of 
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L[i//(t)]  is  shown  to  be  simply  Similarly,  using  the 

convolution  theorem  as  well  as  the  conjugate  symmetry  of 
the  imaginary  part  of  is  seen  to  be  an  inverse 

Hilbert  transform28  of  if/(t ), 

1  f00  if/(r) 

Im[L[^(0]]=--  T—dT=TTlMt)l  (17) 

7T  J  -  oo  t  T 

This  transform,  also  known  as  a  quadrature  filter,  applies  a 
phase  shift  of  7t/2  to  each  frequency  component  of  the  input 
signal. 

Thus,  the  time-domain  reconstruction  formula  can  fi¬ 
nally  be  written 

J  ®(0,a)[ps(0,a,T) 


so  that  f((o)  is  purely  real.  [Similarly,  if  the  incident  pulse 
waveform  is  odd  in  time,  /(<*>)  is  purely  imaginary  and  Eq. 
(18)  can  still  be  employed.] 

However,  supposition  of  a  frequency-independent  phase 
for  f(oy)  does  not  result  in  any  loss  of  generality.  For  any 
linear-phase  signal,  such  that  the  Fourier  transform  has  the 
form 

/(w)  =  l/(w)k'"f>  <o>0,  (21) 

an  additional  delay  term  of  magnitude  £  can  be  applied  to  all 
scattered  signals  to  obtain  the  signals  associated  with  the 
purely-real  spectrum  \f(a>)\.  In  general,  the  scattered  field 
associated  with  a  desired  waveform  f{t)  can  be  determined 
for  an  arbitrary  waveform  u(t)  from  the  deconvolution  op¬ 
eration 


+  /H  1  [p5(  0,  a,  r)]  j  dS  adS  e , 

where 


(18) 


t=R/c0  + 


(a—0)  r 

Co 


(19) 


The  direction-dependent  weight  <E>(0,a),  which  is  the  same 
as  the  “filter”  employed  in  single-frequency  filtered  back- 
propagation,  is  given  for  the  2D  and  3D  cases  by  Eq.  (9). 

Equation  (18)  is  notable  in  several  respects.  First,  it  pro¬ 
vides  a  linearized  reconstruction  that  employs  scattering  in¬ 
formation  from  the  entire  signal  bandwidth  without  any  fre¬ 
quency  decomposition  of  the  scattered  wavefield.  Second, 
the  delay  term  r  corresponds  exactly  to  the  delay  required  to 
construct  a  focus  at  the  point  r  by  delaying  and  summing  the 
scattered  wavefield  ps(0,aj)  for  all  measurement  directions 
0  and  incident-wave  directions  a.  Thus,  the  time-domain 
reconstruction  formula  given  by  Eq.  (18)  can  be  regarded  as 
a  quantitative  generalization  of  confocal  time-domain  syn¬ 
thetic  aperture  imaging,  in  which  signals  are  synthetically 
delayed  and  summed  for  each  transmit/receive  pair  to  focus 
at  the  image  point  of  interest.22  29,30 

A  reconstruction  formula  similar  to,  although  less  gen¬ 
eral  than,  Eq.  (18)  was  independently  derived  in  Ref.  23  for 
the  two-dimensional  inverse  scattering  problem.  In  view  of 
the  present  derivation,  the  method  of  “probing  by  plane 
pulses”  in  Ref.  23  can  be  regarded  to  yield  a  multiple- 
frequency  reconstruction  of  Re[yM(r)],  while  the  present 
method  yields  the  complex  function  yM(  r).  In  Ref.  23,  this 
method  was  proposed  as  a  more  convenient  way  to  imple¬ 
ment  narrow-band  diffraction  tomography.  However,  the  nu¬ 
merical  results  given  below  show  that  the  reconstruction  for¬ 
mula  of  Eq.  (18),  when  directly  implemented  using 
wideband  signals,  provides  considerable  improvement  in  im¬ 
age  quality  over  narrow-band  reconstructions. 

Reconstructions  using  Eq.  (18)  can  be  performed  using 
any  pulse  waveform.  However,  the  frequency  compounding 
defined  by  Eq.  (10)  is  most  straightforwardly  interpreted  if 
the  frequency  weight  g(a>)  has  a  phase  that  is  independent  of 
frequency.  This  criterion  can  be  met,  for  instance,  if  the  in¬ 
cident  pulse  waveform /(r)  is  even  in  time, 


(20) 


[ps(0,a,t)]m= F  1 


/(<«>) 

u((o) 


\ps(0,a,t)]u{  0 


(22) 


For  stable  deconvolution  using  Eq.  (22),  the  desired  /( a>) 
should  not  have  significant  frequency  components  outside 
the  bandwidth  of  u(o)). 


B.  Extensions  to  the  reconstruction  algorithm 

For  large  tissue  structures  at  high  ultrasonic  frequencies, 
weak  scattering  approximations  such  as  the  Bom  approxima¬ 
tion  are  of  limited  validity.  Thus,  for  problems  of  interest  to 
medical  ultrasonic  imaging,  reconstructed  image  quality  can 
be  improved  by  aberration  correction  methods  that  incorpo¬ 
rate  higher-order  scattering  and  propagation  effects.  The 
present  time-domain  reconstruction  formula  (18)  provides  a 
natural  framework  for  quantitative  imaging  with  aberration 
correction.  In  general,  if  the  background  medium  is  known 
or  can  be  estimated,  the  received  scattered  signals  can  be 
processed  to  provide  an  estimate  of  the  scattered  field  that 
would  be  measured  for  the  same  scatterer  within  a  homoge¬ 
neous  background  medium.  This  approach  essentially  re¬ 
moves  higher-order  scattering  effects  from  the  measured  far 
field  scattering,  so  that  a  Bom  inversion  can  be  performed  on 
the  modified  data;  similar  processes  occur  implicitly  in  many 
nonlinear  inverse  scattering  methods.31 

For  example,  a  simple  implementation  of  aberration  cor¬ 
rection  can  be  derived  if  one  makes  the  assumption  that 
background  inhomogeneities  result  only  in  cumulative  de¬ 
lays  (or  advances)  of  the  incident  and  scattered  wavefronts. 
This  crude  model  does  not  include  many  propagation  and 
scattering  effects  important  to  ultrasonic  aberration,  but  has 
been  shown  to  provide  a  reasonable  first  approximation  of 
local  delays  in  wavefronts  propagating  through  large-scale 
tissue  models.32,33  Given  this  approximation,  the  total  delay 
for  an  angle  <fi  and  a  point  position  r  is  given  by 

<M<M  =  fc(^)”1^--,  (23) 

co 

where  the  integral  is  performed  along  the  line  that  joins  the 
spatial  points  r  and  R<f> ,  Aberration-corrected  reconstruc¬ 
tions  can  then  be  performed  using  Eq.  (18)  with  r  replaced 
by  the  corrected  delay  term 
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II.  COMPUTATIONAL  METHODS 


(  or—  0)  •  r 

r— >R/c0~\ - h  <5r(a,r)  +  8t{0,t).  (24) 

co 

Improved  approximations  could  be  obtained  by  application 
of  the  delay  function  8t(<J>, r)  after  numerical  backpropaga- 
tion  of  the  far-field  scattered  wavefronts  through  a  homoge¬ 
neous  medium34,35  or  by  compensation  for  both  delay  and 
amplitude  variations.36,37  More  general,  although  much  more 
computationally  expensive,  aberration  correction  could  also 
be  performed  by  synthetic  focusing  using  full-wave  numeri¬ 
cal  computation  of  acoustic  fields  within  an  estimated  real¬ 
ization  of  the  unknown  medium.  A  method  of  this  kind  has 
been  implemented,  within  the  context  of  a  frequency-domain 
diffraction  tomography  method,  in  Ref.  19. 

The  present  imaging  method  has  been  derived  using 
simplifying  assumptions  including  zero  absorption  and  con¬ 
stant  density  for  the  scattering  medium.  However,  these  as¬ 
sumptions  do  not  substantially  restrict  the  validity  of  the 
method.  For  example,  the  effect  of  absorption  can  be  reduced 
using  time-gain  compensation,  with  or  without  frequency- 
dependent  corrections,38  of  received  scattered  signals  for 
each  transmit/receive  pair.  Such  time-gain  compensation 
could  be  performed  either  using  an  estimated  bulk  attenua¬ 
tion  for  the  medium  (as  with  current  clinical  ultrasound  scan¬ 
ners),  or  by  implementation  of  an  adaptive  attenuation  model 
in  a  manner  similar  to  the  time-shift  compensation  scheme 
discussed  above. 

Inclusion  of  density  variations  as  well  as  sound  speed 
variations  adds  additional  complication  to  the  time-domain 
diffraction  tomography  algorithm  derived  here.  For  single- 
frequency  diffraction  tomography  in  the  presence  of  sound 
speed  and  density  variations,  the  quantity  yB(r,(o)  recon¬ 
structed  by  Eq.  (8)  can  be  shown39  to  provide  an  estimate  of 
a  physical  quantity  that  depends  both  on  sound  speed  varia¬ 
tions  and  density  variations.  In  the  notation  used  here,  this 
quantity  can  be  written 


y'(r)=y(r)  — y(r)yp(r)  +  —  V2yp(r),  (25) 

where  the  density  variation  is  defined  yp-  1  —  p0/p(r). 
Thus,  for  time-domain  reconstructions  of  media  with  density 
variations,  the  reconstruction  formula  of  Eq.  (18)  will  pro¬ 
vide  the  estimate 


y3f(r)asy(r)-y(r)yp(r)  +  -E  V2yp(r),  (26) 

2&o 

where  k0  is  the  wave  number  corresponding  to  the  center 
frequency  of  the  pulse  employed.  For  media  such  as  human 
tissue,  where  density  variations  are  fairly  small  and  abrupt 
density  transitions  are  rare,  the  last  two  terms  of  Eq.  (26)  are 
small  compared  to  y(r),  so  that  the  reconstruction  algorithm 
derived  above  can  still  be  regarded  to  provide  an  image  of 
the  sound-speed  variation  function  y(r).  However,  if  de¬ 
sired,  a  reconstruction  employing  pulses  with  two  distinct 
center  frequencies  could  allow  separation  of  sound  speed  and 
density  variations  by  techniques  similar  to  those  described  in 
Ref.  16  or  39. 


The  time-domain  inverse  scattering  method  described 
above  has  been  tested  with  2D  and  3D  synthetic  data  pre¬ 
pared  using  three  numerical  methods:  a  Bom  approximation 
method  for  point  scatterers  and  3D  slabs,  an  exact  series 
solution  for  cylindrical  inhomogeneities,  and  a  k- space 
method  for  arbitrary  2D  inhomogeneous  media. 

The  time-domain  waveform  employed  for  all  the  com¬ 
putations  reported  here  was 

f(t)  =  cos(  (O0t)  e -  ?2/(2<r2),  (27) 

where  cd0=2  7r/0  for  a  center  frequency  of  /0  and  cr  is  the 
temporal  Gaussian  parameter.  This  waveform  has  the  real, 
even  Fourier  transform 

f((o)=  o)2/2+<T^("+"o>2'2).  (28) 

*  87T 

Values  used  for  the  computations  reported  here  were  /0 
=  2.5  MHz  and  <7=0.25  /zs,  so  that  the  —  6  dB  bandwidth 
of  the  signal  was  1.5  MHz.  These  parameters  correspond 
closely  to  those  of  an  existing  2048-element  ring 
transducer.40 

For  the  case  of  point  scatterers,  the  contrast  function  y 
was  assumed  to  take  the  form 

M 

y(r)  =  2  My^r-rj).  (29) 

Using  the  far-field  form  of  the  2D  Green’s  function  and  ne¬ 
glecting  multiple  scattering,  Eq.  (6)  for  the  scattered  far  field 
can  be  rewritten  as 

P.A  0, a,a))=-k2^j /( w) S  V,ema~e)  'ci  (30) 

for  each  frequency  component  of  interest.  Time-domain 
waveforms  were  synthesized  by  using  Eq.  (30)  for  each  fre¬ 
quency  with  /(<o)>10~3  and  inverting  the  frequency- 
domain  scattered  wavefield  by  a  fast  Fourier  transform  (FFT) 
implementation  of  Eq.  (4).  The  temporal  sampling  rate  em¬ 
ployed  was  10  MHz.  An  analogous  formula,  with  a  different 
multiplicative  constant,  was  also  employed  for  the  3D  case. 

The  Bom  approximation  was  also  used  to  compute 
three-dimensional  scattering  for  slab-shaped  objects  defined 
by  the  equation 

y{r)  =  y0H{ax-\x\)H{ay-\y\)H{az-\z\).  (31) 

For  this  object,  the  linearized  forward  problem  can  be  solved 
analytically.  Under  the  Bom  approximation,  the  frequency- 
domain  scattered  far  field  has  the  form 


Ps(  0,  a,  w)  =  2 /( w)  To  axayazeikR/(  ttR) 


sin [kLx(  a —  6)  •  ej  sin[/:Lv(  a—  0)  •  ev] 
kLx(  a—  0)  •  e*  kLy(  a -  0)  •  ev 


sin [kLz(a—ff)  •  ej 
X  kLz(a—0)-ez 


(32) 


where  ex ,  ev ,  and  ez  represent  unit  vectors  in  the  x ,  y,  and  z 
directions.  The  time  domain  scattered  pressure  ps(0,a,t)  is 
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obtained,  as  for  the  point  scatterer  case  described  above,  by 
inverse  transformation  of  the  frequency-domain  wavefield 
for  all  frequencies  within  the  bandwidth  of  interest. 

For  2D  cylindrical  inhomogeneities,  an  analogous  pro¬ 
cedure  was  followed,  except  that  the  frequency-domain  scat¬ 
tered  wavefield ps(0,a,<o)  was  computed  using  an  exact  se¬ 
ries  solution25  for  each  frequency  component  of  interest.  In 
implementation  of  the  series  solution,  summations  were  trun¬ 
cated  when  the  magnitude  of  a  single  coefficient  dropped 
below  10“ 12  times  the  sum  of  all  coefficients. 

Solutions  were  also  obtained  for  arbitrary  2D  inhomoge¬ 
neous  media  using  a  time-domain  &-space  method.41  Grid 
sizes  of  256X256  points,  a  spatial  step  of  0.0833  mm,  and  a 
time  step  of  0.02734  ps  were  employed.  Scattered  acoustic 
pressure  signals  on  a  circle  of  virtual  receivers  were  recorded 
at  a  sampling  rate  of  9.144  MHz.  The  receiver  circle,  which 
had  a  radius  of  3.0  mm  in  these  computations,  completely 
contained  the  inhomogeneities  used.  Far-field  waveforms 
were  computed  by  Fourier  transforming  the  time-domain 
waveforms  on  the  near-field  measurement  circle,  transform¬ 
ing  these  to  far-field  waveforms  for  each  frequency  using  a 
numerically  exact  transformation  method,16  and  performing 
inverse  Fourier  transformation  to  yield  time-domain  far-field 
waveforms.  All  forward  and  inverse  temporal  Fourier  trans¬ 
forms,  as  well  as  angular  transforms  occurring  in  the  near¬ 
field-far- field  transformation,16  were  performed  by  FFT. 

The  time-domain  imaging  method  was  directly  imple¬ 
mented  using  Eq.  (18),  evaluated  using  straightforward  nu¬ 
merical  integration  over  all  incident-wave  and  measurement 
directions  employed.  The  reconstruction  formula  employed 
can  be  explicitly  written  as 

1  C27T  C  2tt  j 

|sin(a—  0)|l  ps{9,a,r) 

A  2d  Jo  Jo  \ 


FIG.  2.  Point-spread  functions  for  time-domain  and  single -frequency  dif¬ 
fraction  tomography  methods.  In  each  panel,  the  vertical  scale  corresponds 
to  the  relative  amplitude  of  the  reconstructed  contrast  y(r),  while  the  hori¬ 
zontal  scale  corresponds  to  number  of  wavelengths  at  the  center  frequency, 
(a)  Two-dimensional  case,  (b)  Three-dimensional  case. 
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c0 

for  the  2D  case,  where  a  and  6  are  the  angles  corresponding 
to  the  direction  vectors  a  and  0 ,  and  as 

1  C2Tt  f  77"  r  2ir  C  TT  I 

7m(  r)=77—  \a-mps{e,a,T) 
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for  the  3D  case,  where  0a  and  a  are  direction  angles  for 
the  incident- wave  direction  a  and  0  e  and  <J>  0  are  direction 


angles  for  the  measurement  direction  0 .  For  each  case,  the 
normalization  factor  N  was  determined  from  Eq.  (11)  with 
g{o))—f{o))l &{(*))  and  /x(o>)  given  by  Eq.  (9).  Before 
evaluation  of  the  argument  r  for  each  signal,  the  time- 
domain  waveforms  were  resampled  at  a  sampling  rate  of  16 
times  the  original  rate.  This  resampling  was  performed  using 
FFT-based  Fourier  interpolation.  The  inverse  Hilbert  trans¬ 
form  was  performed  for  each  signal  using  an  FFT  implemen¬ 
tation  of  Eq.  (16).  Values  of  the  pressure  signals  at  the  time 
r  were  then  determined  using  linear  interpolation  between 
samples  of  the  resampled  waveforms.  The  integrals  of  Eqs. 
(33)  and  (34)  were  implemented  using  discrete  summation 
over  all  transmission  and  measurement  directions  employed. 

Computations  were  also  performed  using  the  time- 
domain  diffraction  tomography  algorithm  for  limited- 
aperture  data.  For  these  reconstructions,  the  integrals  of  Eq. 
(33)  were  evaluated  only  for  angles  corresponding  to  trans¬ 
mitters  and  receivers  within  a  specified  aperture  of  angular 
width  <£ap,  i.e., 

|a|^ap/2,  1 0—  tt\  0ap/2.  (35) 

Use  of  a  small  value  for  <£ap  corresponds  to  use  of  a  small 
aperture  in  pulse-echo  mode. 
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HI.  NUMERICAL  RESULTS 


Two-dimensional  and  three-dimensional  point-spread 
functions  (PSF)  for  the  present  time-domain  diffraction  to¬ 
mography  method  are  illustrated  in  Fig.  2.  The  time-domain 
reconstructions  shown  here,  like  the  other  time-domain  re¬ 
constructions  shown  in  this  paper,  were  obtained  using  a 
incident  pulse  of  center  frequency  2.5  MHz  and  a  Gaussian 
envelope  corresponding  to  a  —  6  dB  bandwidth  of  1 .5  MHz. 
Point-spread  functions  were  determined  by  reconstructing  a 
point  scatterer  located  at  the  origin.  For  the  2D  case,  in 
which  the  point  scatterer  can  be  regarded  as  a  thin  wire, 
synthetic  scattering  data  was  obtained  using  the  Bom  ap¬ 
proximation  method  outlined  above  for  16  incident- wave  di¬ 
rections  and  64  measurement  directions.  The  3D  time- 
domain  reconstruction  was  obtained  using  Bom  data  for  72 
incident-wave  directions  and  288  measurement  directions, 
each  evenly  spaced  on  a  rectangular  grid  defined  by  the 
angles  ©  and  <t>.  For  comparison,  analogous  point-spread 
functions  are  also  shown  for  standard  frequency-domain  dif¬ 
fraction  tomography  reconstructions  using  single-frequency 
(2.5  MHz)  data. 

For  the  2D  case  illustrated  in  Fig.  2,  the  time-domain 
PSF  has  a  slightly  narrower  peak,  indicating  that  point  reso¬ 
lution  has  been  slightly  improved  by  the  increased  band¬ 
width  employed  in  the  time  domain  method.  More  signifi¬ 
cantly,  sidelobes  of  the  time-domain  PSF  are  significantly 
smaller  than  those  for  the  single-frequency  PSF  (the  first 
sidelobe  is  reduced  by  7  dB,  while  the  second  is  reduced  by 
19  dB),  so  that  contrast  resolution  for  time-domain  diffrac¬ 
tion  tomography  is  seen  to  be  much  higher  than  for  single- 
frequency  diffraction  tomography.  For  the  3D  case,  the  time- 
domain  reconstruction  shows  a  much  more  dramatic 
improvement  over  the  single-frequency  reconstruction.  In 
this  case,  the  time-domain  solution  shows  significant  in¬ 
creases  in  both  the  point  resolution  (PSF  width  at  half¬ 
maximum  reduced  by  27%)  and  contrast  resolution  (first 
sidelobe  reduced  by  13  dB  and  second  sidelobe  reduced  by 
18  dB).  Furthermore,  a  comparison  of  the  PSFs  for  2D  and 
3D  time-domain  reconstruction  indicates  that  much  higher 
image  quality  is  achievable  for  3D  time-domain  imaging 
than  for  the  2D  case.  This  increase  in  image  quality  suggests 
that  the  time-domain  diffraction  tomography  method  pro¬ 
posed  here  may  benefit  from  the  overdetermined  nature  of 
the  general  wideband  3D  inverse  scattering  problem.42,43 

The  effect  of  transmit  and  receive  aperture  characteris¬ 
tics  on  image  quality  is  illustrated  in  Fig.  3.  Panels  (a)  and 
(b)  of  Fig.  3  show  the  point-spread  function  for  a  number  of 
aperture  configurations,  each  employing  64  measurement  di¬ 
rections.  Figure  3(a)  shows  the  point-spread  function  for  re¬ 
constructions  obtained  using  1,  4,  8,  and  16  incident- wave 
directions.  The  point  scatterer  is  clearly  imaged  even  for  the 
reconstruction  using  one  incident-wave  direction.  Optimal 
image  quality  (indistinguishable  from  reconstructions  with 
64  incident- wave  directions)  is  obtained  for  16  incident- 
wave  directions,  so  that  scattering  data  obtained  using  one 
incident-wave  direction  for  each  group  of  four  measurement 
directions  appears  to  be  sufficient  for  the  present  reconstruc¬ 
tion  method. 

The  effect  of  limited  view  range  on  the  point  spread 
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FIG.  3.  Effect  of  aperture  characteristics  on  image  quality.  Each  panel 
shows  the  real  part  of  a  time-domain  reconstruction.  Re  [yM],  on  a  linear 
grayscale  with  white  representing  the  maximum  amplitude  of  |yw(r)|  and 
black  represents  —  1  times  the  maximum  amplitude,  (a)  Point-spread  func¬ 
tions  for  the  same  waveform  parameters  as  Fig.  2.  Each  panel  shows  an  area 
of  0.6 X  0.6  mm2,  corresponding  to  one  square  wavelength  at  the  center 
frequency.  Left  to  right:  1,4,  8,  and  16  incident- wave  directions,  (b)  Point- 
spread  functions  for  aperture  sizes  of  7t/2,  7 r,  3  it/2,  and  2 1 r  radians,  format 
as  in  previous  panel,  (c)  Real  parts  of  reconstructions  for  a  homogeneous 
cylinder  {a—  1.0  mm,  y=0.02).  The  area  shown  in  each  panel  is  2.0X2.0 
mm2. Left  to  right:  aperture  sizes  of  7t/2,  tt,  3it/2 ,  and  27t  radians. 

function  is  also  illustrated  in  Fig.  3.  Panel  (b)  shows  the 
point-spread  function  for  four  differently  limited  apertures, 
while  panel  (c)  shows  reconstructions  of  a  homogeneous  cyl¬ 
inder  (a=  1.0  mm,  y=0.02)  for  the  same  apertures.  In  each 
case,  limitation  of  the  transmit  and  receive  apertures  to 
angles  near  the  backscatter  direction  (aperture  size  it  IT)  re¬ 
sults  in  images  that  resemble  a  conventional  B-scans.  Use  of 
apertures  corresponding  to  pulse-echo  mode  in  the  large- 
aperture  limit  (aperture  size  i t)  yield  higher  resolution  in  all 
directions.  Using  three-fourths  of  a  circular  aperture  (size 
3  7t/2)  yields  image  quality  close  to  that  for  the  full  aperture 
(27r)  case.  The  characteristics  of  all  these  images  result  from 
the  set  of  spatial-frequency  vectors  interrogated  by  each 
group  of  scattering  measurements.1  Apertures  with  only  a 
limited  range  of  transmit  and  receive  directions  [e.g.,  the 
“b-scan”  apertures  shown  in  the  first  column  of  panels  (b) 
and  (c)]  provide  only  information  corresponding  to  large 
spatial  frequency  vectors  oriented  nearly  on-axis,  so  that 
such  images  mainly  show  those  edges  that  are  nearly  perpen¬ 
dicular  to  the  axis  of  the  aperture. 

Reconstructions  performed  using  exact  solutions  for 
scattering  from  cylindrical  inhomogeneities  provide  a 
straightforward  means  to  assess  the  accuracy  of  the  time- 
domain  scattering  method  for  a  range  of  object  sizes  and 
contrasts.  A  number  of  example  reconstructions  are  shown  in 
Figs.  4  and  5.  The  number  of  measurement  directions  for  all 
cylinder  reconstructions  was  chosen  based  on  an  empirical 
test  of  the  number  required  for  a  satisfactory  image  of  a 
homogeneous  cylinder;  for  a  cylinder  of  radius  1  mm,  the 
required  number  of  measurement  directions  was  determined 
to  be  approximately  96.  Based  on  spatial-frequency  sampling 
considerations,  the  number  of  measurement  directions  was 
increased  in  proportion  to  the  size  of  the  inhomogeneous 
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FIG.  4.  Cross  sections  of  reconstructed  contrast  functions  y(r)  for  a  cylinder  of  radius  1  mm,  using  time-domain  (TD)  and  single-frequency  (SF)  diffraction 
tomography.  Waveform  parameters  are  as  in  Fig.  1.  (a)  y=0.02.  (b)  y=0.04.  (c)  y— 0.06.  (d)  y=0.08. 
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region  to  be  reconstructed.  Since  the  results  shown  in  Fig.  3 
indicate  that  considerably  fewer  incident-wave  directions 
than  measurement  directions  are  needed,  the  number  of  inci¬ 
dent  directions  was  chosen  to  be  one-quarter  the  number  of 
measurement  directions  in  each  case. 

Cross  sections  of  time-domain  and  single-frequency  re¬ 
constructions,  plotted  in  Fig.  4,  show  the  relative  accuracy  of 
each  reconstruction  method  for  a  cylinder  of  1-mm  radius 
and  purely  real  contrast  ranging  from  y—  0.02  to  y=0.08. 
For  the  synthetic  scattering  data  in  each  case,  96  measure¬ 
ment  directions  and  24  incident-wave  directions  were  em¬ 
ployed.  The  time-domain  reconstructions  show  improvement 
over  the  single-frequency  reconstructions  both  in  improved 
contrast  resolution  (smaller  sidelobes  outside  the  support  of 
the  cylinder)  and  in  decreased  ringing  (Gibbs  phenomenon) 
artifacts  within  the  support  of  the  cylinder.  However,  for 
increasing  contrast  values,  both  methods  show  similar  in¬ 
creases  in  phase  error,  as  indicated  by  increased  imaginary 
parts  of  the  reconstructed  contrast.  This  error  results  from  the 
Bom  approximation,  which  is  based  on  the  assumption  that 
the  incident  wave  propagates  through  the  inhomogeneous 
medium  without  distortion.  Perturbations  in  the  local  arrival 
time  of  the  incident  wavefront,  which  are  more  severe  for 
higher  contrasts  and  larger  inhomogeneities,  can  result  in  a 
scattered  field  that  is  phase  shifted  relative  to  the  ideal  case 
assumed  in  the  Born  approximation;  linear  inversion  of  this 
phase-distorted  data  naturally  results  in  a  phase-distorted  re¬ 
construction  of  the  scattering  medium.  (A  complementary 
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explanation  of  this  phase  error,  based  on  the  unitarity  of  the 
scattering  operator,  is  given  in  Ref.  19.) 

A  test  of  image  fidelity  for  the  time-domain  reconstruc¬ 
tion  method  is  shown  in  Fig.  5.  The  real  parts  of  time- 
domain  reconstructions  are  shown  as  grayscale  images  for 
homogeneous  cylinders  with  radii  between  1  and  4  mm  and 
contrasts  between  y=0.02  and  y=0.08.  The  number  of 
measurement  directions  employed  for  the  synthetic  scatter¬ 
ing  data  was  96  for  the  1-mm  radius  cylinders,  192  for  the 
2-mm  cylinders,  288  for  the  3-mm  cylinders,  and  384  for  the 
4-mm  cylinders.  In  each  case,  four  incident-wave  directions 
per  measurement  direction  were  used.  The  first  row  of  this 
figure  corresponds  to  the  time-domain  reconstructions  shown 
in  Fig.  4. 

The  images  shown  in  Fig.  5  provide  a  basis  for  evaluat¬ 
ing  the  ability  of  the  present  time-domain  diffraction  tomog¬ 
raphy  method  to  image  homogeneous  objects  of  various 
sizes  and  contrasts.  In  this  figure,  images  of  Re[yM]  show 
uniform  quality  for  small  cylinder  sizes  and  contrasts,  but 
poorer  image  quality  for  larger  sizes  and  contrasts.  For  the 
largest  size  and  contrast  employed  (a  —  4.0  mm,  y=0.08), 
the  reconstruction  primarily  shows  the  edges  of  the  cylinder 
and  fails  to  image  the  interior.  Particularly  notable  is  that  the 
“matrix”  of  images  in  Fig.  5  is  nearly  diagonal;  that  is,  a 
linear  increase  in  object  contrast  causes  image  degradation 
comparable  to  a  corresponding  linear  increase  in  object  size. 
Thus,  a  nondimensional  parameter  directly  relevant  to  image 
quality  for  homogeneous  objects  is  ka  y,  where  k  is  a  domi- 
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FIG.  5.  Images  of  time-domain  reconstructions  for  cylinders  of  varying 
radius  a  and  contrast  y.  Each  panel  shows  the  real  part  of  the  reconstructed 
contrast,  Re[y^(r)],  for  a  pulse  of  center  frequency  2.5  MHz  and  -6  dB 
bandwidth  1.5  MHz.  The  area  shown  in  each  panel  is  2a  X  2a.  All  images 
are  shown  on  a  linear,  bipolar  gray  scale  where  white  represents  the  maxi¬ 
mum  amplitude  of  |yM(r)|  and  black  represents  —1  times  the  maximum 
amplitude. 

nant  wave  number,  a  is  the  object  radius,  and  y  is  the  object 
contrast.  Using  the  wave  number  k0=  10.472  rad/mm  corre¬ 
sponding  to  the  center  frequency  of  2.5  MHz  and  a  sound 
speed  of  L5  mm/ju,s,  the  reconstructions  shown  in  Fig.  5 
indicate  that  the  interior  of  the  cylinder  is  imaged  satisfacto¬ 
rily  for  the  approximate  range  ka  y< 2.5.  This  result  is  con¬ 
sistent  with  a  previous  study  of  single-frequency  diffraction 
tomography,  in  which  adequate  Bom  reconstructions  of  cyl¬ 
inders  were  obtained  for  the  parameter  range  ka  y^  2.2. 44 

Reconstructions  for  several  scattering  objects  without 
special  symmetry  are  shown  in  Fig.  6.  All  of  these  recon¬ 
structions  were  performed  using  synthetic  data  produced  by 
the  fc-space  method  described  in  Ref.  41.  Synthetic  scattering 
data  were  computed  for  64  incident-wave  directions  and  256 
measurement  directions  in  each  case.  The  first  panel  shows  a 
reconstruction  of  a  cylinder  of  radius  2.5  mm  and  contrast 
y  =  -0.0295  with  an  internal  cylinder  of  radius  0.2  mm  and 
contrast  y—  0.0632.  These  contrast  values  correspond,  based 
on  tissue  parameters  given  in  Ref.  32,  to  the  sound-speed 
contrasts  of  human  skeletal  muscle  for  the  outer  cylinder  and 
of  human  fat  for  the  inner  cylinder.  The  second  panel  shows 
a  reconstruction  of  a  2.5-mm-radius  cylinder  with  random 
internal  structure.  The  third  reconstruction  shown  employed 
a  portion  of  a  chest  wall  tissue  map  from  Ref.  45.  In  this 
case,  the  synthetic  data  was  obtained  using  a  tissue  model45 
that  incorporates  both  sound  speed  and  density  variations,  so 
that  the  reconstructed  quantity  is  given  by  Eq.  (26).  In  Fig. 
6(c),  black  denotes  connective  tissue  (y=  —  0.1073,  yp 
-0.1134),  dark  gray  denotes  muscle  (y~  -0.0295,  yp 
=  0.0543),  and  light  gray  denotes  fat  (y=0.0632,  yp 
=  —0.0453). 

The  real  part  of  each  reconstruction  in  Fig.  6  shows 
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FIG.  6.  Time-domain  reconstructions  from  full-wave  synthetic  data  for  three 
arbitrary  scattering  objects.  Each  row  shows  the  actual  (purely  real)  contrast 
function  y  together  with  the  real  and  imaginary  parts  of  the  reconstructed 
contrast  function  yM ,  using  the  same  linear  bipolar  gray  scale  for  each 
panel.  Each  panel  shows  a  reconstruction  area  of  5x5  mm2,  (a)  Cylinder, 
radius  2.5  mm,  with  an  internal  cylinder  of  radius  0.2  mm.  (b)  Cylinder, 
radius  2.5  mm,  with  random  internal  structure,  (c)  Tissue  structure,  with 
variable  sound  speed  and  density,  from  chest  wall  cross  section  5L  in 
Ref.  45. 


good  image  quality,  with  high  resolution  and  very  little  evi¬ 
dence  of  artifacts.  Particularly  notable  is  the  accurately  de¬ 
tailed  imaging  of  internal  structure  for  the  random  cylinder 
and  the  chest  wall  cross  section.  As  expected,  the  density 
variations  present  in  the  chest  wall  cross  section  have  not 
greatly  affected  the  image  appearance;  there  is,  however,  a 
slight  edge  enhancement,  associated  with  the  Laplacian  term 
in  Eq.  (26),  at  boundaries  between  tissue  regions.  Also  no¬ 
table  is  the  nearly  complete  absence  of  any  artifacts  outside 
the  scatterer  in  each  case;  this  result  indicates  that  high  con¬ 
trast  resolution  has  been  achieved.  However,  in  each  case, 
the  imaginary  part  of  the  reconstruction  is  nonzero,  indicat¬ 
ing  that  the  Bom  approximation  is  not  fully  applicable.  The 
imaginary  parts  of  each  reconstruction  are,  however,  small 
compared  to  the  real  parts.  Thus,  simple  aberration  correc¬ 
tion  methods  [of  which  one  example  is  given  by  Eq.  (24)] 
could  substantially  reduce  this  phase  error,  as  for  multiple- 
frequency  diffraction  tomography  in  Ref.  19. 

Three-dimensional  reconstructions  of  a  homogeneous 
slab  are  shown  in  Fig.  7.  The  scatterer  is  characterized  by  Eq. 
(31)  with  y0  =  0.01,  a*  =  0.5  mm,  ay=l.O  mm,  and  az=  1.5 
mm.  Synthetic  data  was  computed  using  Eq.  (34)  for  288 
incident- wave  directions  and  1152  measurement  directions, 
each  evenly  spaced  in  the  angles  <f>  and  ©.  Signal  param¬ 
eters  were  as  for  the  examples  above,  except  that  the  initial 
sampling  rate  for  the  time-domain  signals  was  9.0  MHz.  Iso¬ 
surface  renderings  of  the  real  part  of  yM  are  shown  for  the 
surfaces  yM  =  0.0025.  Since  the  scattering  data  were  ob¬ 
tained  using  a  Born  approximation  for  the  3D  case,  the 
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FIG.  7.  Three-dimensional  reconstructions  of  a  uniform  slab  with  contrast 
y=0.01.  Each  reconstruction  shows  an  isosurface  rendering  of  the  surface 
yM= 0.0025.  Left:  single-frequency  reconstruction.  Right:  time-domain  re¬ 
construction. 

imaginary  part  of  each  reconstruction  is  identically  zero  for 
both  reconstructions.  Consistent  with  the  point-spread  func¬ 
tions  shown  in  Fig.  2,  the  time-domain  reconstruction  is 
much  more  accurate  than  the  single-frequency  reconstruc¬ 
tion.  While  the  single-frequency  reconstruction  shows  an  er¬ 
roneously  rippled  surface,  the  time-domain  reconstruction  is 
smooth.  The  time-domain  reconstruction  is  nearly  identical 
to  the  original  object  except  for  some  rounding  of  the  sharp 
edges  due  to  the  limited  high-frequency  content  of  the  signal 
employed.  The  length  scale  of  the  rounded  edges  is  on  the 
order  of  one-half  the  wavelength  of  the  highest  frequency  in 
the  pulse,  i.e.,  about  0.2  mm  for  the  —  6-dB  cutoff  of  3.25 
MHz. 

Since  three-dimensional  inverse  scattering  is  a  computa¬ 
tionally  demanding  problem,  comparison  of  computational 
efficiency  for  single-frequency  and  time-domain  methods  is 
of  interest.  For  both  reconstructions  shown  in  Fig.  7,  identi¬ 
cal  discretizations  of  the  reconstructed  medium  were  em¬ 
ployed.  Both  computations  included  solution  of  the  appli¬ 
cable  linearized  forward  problem  as  well  as  the  inverse 
problem.  Nonetheless,  the  time-domain  method  was  more 
efficient  than  the  single-frequency  method;  the  total  CPU 
time  required  on  a  200-MHz  AMD  K6  processor  was  133.3 
CPU  min  for  the  time-domain  method  and  287.4  CPU  min 
for  the  single- frequency  method.  This  gain  in  efficiency  was 
possible  because  the  greatest  computational  expense  oc¬ 
curred  in  the  “backpropagation”  of  the  signals  for  each  re¬ 
construction  point.  For  the  single-frequency  method,  this 
step  required  evaluation  of  complex  exponentials  for  each 
incident-wave  direction,  measurement  direction,  and  spatial 
point.  For  the  time-domain  method,  however,  the  computa¬ 
tionally  intensive  steps  (including  the  forward  problem  solu¬ 
tion  and  Fourier  interpolation  of  the  scattered  signals) 
needed  only  to  be  performed  once  for  each  transmit/receive 
pair.  For  the  backpropagation  step,  performed  at  each  point 
in  the  3D  spatial  grid,  the  time-domain  reconstruction 
method  required  only  linear  interpolation  of  the  oversampled 
farfield  pressure  waveforms. 

IV.  CONCLUSIONS 

A  new  method  for  time-domain  ultrasound  diffraction 
tomography  has  been  presented.  The  method  provides  quan¬ 
titative  images  of  sound  speed  variations  in  unknown  media; 
when  two  pulse  center  frequencies  are  employed,  the  method 
is  also  capable  of  imaging  density  variations.  Reconstruc¬ 
tions  performed  using  this  method  are  equivalent  to  multiple- 
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frequency  reconstructions  using  filtered  backpropagation,  buf 
can  be  obtained  with  much  greater  efficiency. 

The  time-domain  reconstruction  algorithm  has  been  de¬ 
rived  as  a  simple  filtered  delay -and-sum  operation  applied  to 
far-field  scattered  signals.  This  algorithm  is  closely  related  to 
time-domain  confocal  synthetic  aperture  imaging,  so  that  it 
can  be  considered  a  generalization  of  imaging  algorithms 
employed  in  current  clinical  instruments.  The  simplicity  of 
the  imaging  algorithm  allows  straightforward  addition  of  fea¬ 
tures  such  as  time-gain  compensation  and  aberration  correc¬ 
tion. 

Numerical  results  obtained  using  synthetic  data  for  2D 
and  3D  scattering  objects  show  that  the  time-domain  method 
can  yield  significantly  higher  image  quality  (and,  in  some 
cases,  also  greater  efficiency)  than  single-frequency  diffrac¬ 
tion  tomography.  Quantitative  reconstructions,  obtained  us¬ 
ing  signal  parameters  comparable  to  those  for  present-day 
clinical  instruments,  show  accurate  imaging  of  objects  with 
simple  deterministic  structure,  random  internal  structure,  and 
structure  based  on  a  cross-sectional  tissue  model.  The 
method  is  hoped  to  be  useful  for  diagnostic  imaging  prob¬ 
lems  such  as  the  detection  and  characterization  of  lesions  in 
ultrasonic  mammography. 
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Abstract — A  quantitative  ultrasonic  imaging  method  employing 
time-domain  scattering  data  is  presented.  This  method  provides  to¬ 
mographic  images  of  medium  properties  such  as  the  sound  speed 
contrast;  these  images  are  equivalent  to  multiple-frequency  filtered- 
backpropagation  reconstructions  using  all  frequencies  within  the 
bandwidth  of  the  incident  pulse  employed.  However,  image  synthesis  is 
performed  directly  in  the  time  domain  using  coherent  combination  of 
farfield  scattered  pressure  waveforms,  delayed  and  summed  to  numeri¬ 
cally  focus  on  the  unknown  medium.  The  time-domain  method  is  more 
efficient  than  multiple-frequency  diffraction  tomography  methods,  and 
can,  in  some  cases,  be  more  efficient  than  single-frequency  diffrac¬ 
tion  tomography.  Example  reconstructions,  obtained  using  synthetic 
data  for  two-dimensional  and  three-dimensional  scattering  of  wide¬ 
band  pulses  as  well  as  measured  scattering  data  from  a  2048-element 
ring  transducer,  show  that  the  time-domain  reconstruction  method 
provides  image  quality  superior  to  single-frequency  reconstructions  for 
objects  of  size  and  contrast  relevant  to  medical  imaging  problems  such 
as  ultrasonic  mammography.  The  present  method  is  closely  related  to 
existing  synthetic-aperture  imaging  methods  such  as  those  employed 
in  clinical  ultrasound  scanners.  Thus,  the  new  method  can  be  extended 
to  incorporate  available  image-enhancement  techniques  such  as  time- 
gain  compensation  to  correct  for  medium  absorption  and  aberration 
correction  methods  to  reduce  error  associated  with  weak  scattering  ap¬ 
proximations. 

I.  Introduction 

Quantitative  imaging  of  tissue  properties  is  a  potentially 
useful  technique  for  diagnosis  of  cancer  and  other  disease. 
Inverse  scattering  methods  such  as  diffraction  tomography 
can  provide  quantitative  reconstruction  of  tissue  properties 
including  sound  speed,  density,  and  absorption.  However, 
although  previous  inverse  scattering  methods  have  achieved 
high  resolution  and  quantitative  accuracy,  such  methods 
have  not  yet  been  incorporated  into  commercially  successful 
medical  ultrasound  imaging  systems.  Previous  methods  of 
diffraction  tomography  have  usually  been  based  on  single¬ 
frequency  scattering,  while  current  diagnostic  ultrasound 
scanners  employ  wideband  time-domain  signals.  The  use 
of  wideband  information  in  image  reconstruction  is  known 
to  provide  increased  point  and  contrast  resolution,  both  of 
which  are  important  for  medical  diagnosis. 

Relatively  few  previous  workers  have  investigated  direct 
use  of  wideband  scattering  data  for  inverse  scattering  meth¬ 
ods  analogous  to  single-frequency  diffraction  tomography. 
A  review  of  several  approaches  is  given  in  Ref.  [1],  includ¬ 
ing  linear  and  nonlinear  diffraction  tomography  methods  us¬ 
ing  scattering  data  for  a  number  of  discrete  frequencies  [2]— 
[4],  a  direct  (but  not  completely  general)  time-domain  re¬ 
construction  algorithm  [5],  and  an  extension  of  the  eigen¬ 
function  method  from  Ref.  [4]  to  use  the  full  bandwidth  of 
the  incident  pulse  waveform  [6]. 

Recently,  a  new  approach  to  wideband  quantitative  imag¬ 
ing  has  been  offered:  a  time-domain  inverse  scattering 


method  that  overcomes  some  of  the  limitations  of  previ¬ 
ous  frequency-domain  and  time-domain  quantitative  imag¬ 
ing  methods  [1].  In  this  paper,  the  new  time-domain  diffrac¬ 
tion  tomography  algorithm  is  briefly  reviewed.  The  capabil¬ 
ities  of  the  method  are  demonstrated  using  simulated  recon¬ 
structions  of  two-dimensional  and  three-dimensional  scat¬ 
tered.  The  practical  capability  of  the  method  for  ultra¬ 
sonic  mammography  is  then  illustrated  by  reconstructions  of 
tissue-mimicking  phantoms  from  scattering  data  measured 
by  a  2.5  MHz,  2048-element  ring  transducer. 

II.  Theory 

A  new  time-domain  inverse  scattering  algorithm,  appli¬ 
cable  to  quantitative  imaging  of  tissue  and  other  inhomoge¬ 
neous  media,  is  derived  in  Ref.  [1]  and  summarized  briefly 
below.  The  medium  is  modeled  as  a  fluid  medium  defined 
by  the  sound  speed  contrast  function  7(1-)  =  <^/c( r)2  -  1, 
where  co  is  a  background  sound  speed  and  c(r)  is  the 
spatially-dependent  sound  speed  defined  at  all  points  r.  For 
the  scope  of  the  present  paper,  the  medium  is  assumed  to 
have  constant  density,  no  absorption,  and  weak  scattering 
characteristics;  extensions  to  the  reconstruction  algorithm 
that  overcome  these  limiting  assumptions  are  discussed  in 
Ref.  [1]. 

The  medium  is  subjected  to  a  pulsatile  plane  wave  of  the 
form  Pmc(r,<M)  =  fit  -  r  •  a/co),  where  a  is  a  unit 
vector  in  the  direction  of  propagation,  /  is  the  time-domain 
waveform,  and  co  is  the  background  sound  speed.  The  scat¬ 
tered  wavefieldp*  (0,  a,  t)  is  measured  at  a  fixed  radius  R  in 
the  farfield,  where  6  corresponds  to  the  direction  unit  vector 
of  a  receiving  transducer  element.  (Alternatively,  if  scat¬ 
tering  measurements  are  made  in  the  nearfield,  the  farfield 
acoustic  pressure  can  be  computed  using  exact  transforms 
that  represent  propagation  through  a  homogeneous  medium 
[2].)  The  farfield  scattered  pressure,  when  specified  for  all 
incident-wave  directions  a,  measurement  directions  6 ,  and 
times  t ,  comprises  the  data  set  to  be  used  for  reconstruction 
of  the  unknown  medium.  The  inverse  scattering  problem  is 
to  reconstruct  the  unknown  medium  contrast  7(1*)  using  the 
scattered  field  p8  ( 0 ,  a,  a;)  measured  at  a  fixed  radius  R. 

The  starting  point  for  the  present  time-domain  inverse 
scattering  method  is  single-frequency  filtered  backpropaga- 
tion  [2],  [7],  [8].  Under  the  assumption  of  weak  scattering, 
such  that  the  Bom  approximation  holds,  the  solution  to  the 
single-frequency  inverse  scattering  problem  is  given  by  the 
formula 

7B(r,w)  =  —  ff  $(0,a)p,{0,a,w) 

f(u)  JJ 

x  eik(9~a]  r dSadS(,,  where  (1) 


1999  IEEE  Ultrasonics  Symposium 


(c)  1999  IEEE 


0-7803-5725-6/99/$1 0.00 


A (w)  =  Y  $(0,a)  =  I  sin(0  -  a)|  in  2D,  and 
AM  =  ||,  •(*,<*)  =  |0-or|  in  3D.  (2) 

Each  surface  integral  in  Eq.  (1)  is  performed  over  the  en¬ 
tire  measurement  circle  for  the  2D  case  and  over  the  entire 
measurement  sphere  for  the  3D  case.  Equation  (1)  provides 
an  exact  solution  to  the  linearized  inverse  scattering  prob¬ 
lem  for  a  single  frequency  component  of  the  scattered  wave- 
field  p6  ( 0 ,  a,  t) .  The  resulting  reconstruction,  jB  (r,  u>),  has 
spatial  frequency  content  limited  by  the  “Ewald  sphere”  of 
radius  2k  in  wavespace  [9], 

To  improve  upon  the  single-frequency  formulas  specified 
by  Eq.  (1),  one  can  extend  the  spatial-frequency  content  of 
reconstructions  by  exploiting  wideband  scattering  informa¬ 
tion.  The  method  outlined  here  synthesizes  a  “multiple- 
frequency”  reconstruction  7 m(f)  by  formally  integrating 
single- frequency  reconstructions  7 B  (r ,  oj)  over  a  range  of 
frequencies  oj.  A  general  formula  for  this  approach  is 


,  _  Jn°°  ffM  Jb  (t,  to)  dip 

T)~  fo°°sMdu 


where  g(oj)  is  an  appropriate  frequency-dependent  weight¬ 
ing  function.  In  practice,  the  weighting  function  g(oj)  is  cho¬ 
sen  to  be  bandlimited  because  (for  a  given  set  of  physical 
scattering  measurements)  the  frequency-dependent  contrast 
7b  (r,  oj)  can  only  be  reliably  reconstructed  for  a  finite  range 
of  frequencies  oj  associated  with  the  spectra  of  the  incident 
waves  employed.  Thus,  the  integrands  in  Eq.  (3)  are  nonzero 
only  over  the  support  of  g(co)  and  the  corresponding  inte¬ 
grals  are  finite. 

If  the  frequency  weighting  function  is  now  specified 
to  incorporate  the  incident-pulse  spectrum  as  well  as  the 
frequency-  and  dimension- dependent  coefficient  fi{oj),  such 
that  g(uj)  =  /M/£(u),  Eq.  (3)  reduces  to  the  form  [1] 

7Af(r)  =  -t  JJ  $(0,a)  (p$(0,a,r) 

+  i  H-1  [p, (0, a,  t)]  )  dSa  dSg ,  where  (4) 

r  =  fi/co  +  ^  ~  — — ,  N  =  2  [  g(u)<L>, 

Co  Jo 

and  H"1  is  the  inverse  Hilbert  transform,  also  known  as  a 
quadrature  filter. 

Equation  (4)  is  notable  in  several  respects.  First,  it  pro¬ 
vides  a  linearized  reconstruction  that  employs  scattering  in¬ 
formation  from  the  entire  signal  bandwidth  without  any  fre¬ 
quency  decomposition  of  the  scattered  wavefield.  Second, 
the  delay  term  r  corresponds  exactly  to  the  delay  required 
to  construct  a  focus  at  the  point  r  by  delaying  and  sum¬ 
ming  the  scattered  wavefield  p6(0,ayt)  for  all  measure¬ 
ment  directions  6  and  incident- wave  directions  a.  Thus,  the 
time-domain  reconstruction  formula  given  by  Eq.  (4)  can  be 
regarded  as  a  quantitative  generalization  of  confocal  time- 
domain  synthetic  aperture  imaging  (e.g.y  the  “gold  standard” 
beamformerof  Ref.  [10]),  in  which  signals  are  synthetically 
delayed  and  summed  for  each  transmit/receive  pair  to  focus 
at  the  image  point  of  interest. 


Fig.  1.  Point-spread  function  for  three-dimensional  time-domain  and 
single-frequency  diffraction  tomography  methods.  The  vertical  scale  cor¬ 
responds  to  the  relative  amplitude  of  the  reconstructed  contrast  7(1*),  while 
the  horizontal  scale  corresponds  to  number  of  wavelengths  at  the  center 
frequency. 


ID.  Simulations 

Below,  the  time-domain  diffraction  tomography  method 
of  Ref.  [1]  is  illustrated  using  results  of  simulation  tests 
with  2D  and  3D  synthetic  data.  The  synthetic  scattering 
data  employed  were  obtained  using  a  Bom  approximation 
method  for  point  scatterers  and  3D  slabs,  and  a  fc-space 
method  [11]  for  arbitrary  2D  inhomogeneous  media.  Ad¬ 
ditional  results,  presented  in  Ref.  [1],  show  reconstructions 
performed  using  using  exact  time-domain  solutions  for  scat¬ 
tering  from  compressible  cylinders  as  well  as  reconstruc¬ 
tions  from  limited-aperture  data.  The  time-domain  wave¬ 
form  employed  for  all  the  simulations  reported  here  was 

f(t)  =  cos(27r/o£)  e_*2/(2cr2\  with  fo  =  2.5  MHz  and 
cr  =  0.25  ^s,  so  tnat  the  -6  dB  bandwidth  of  the  signal 
was  1 .5  MHz.  These  parameters  correspond  closely  to  those 
of  the  ring  transducer  used  in  the  measurements  reported  in 
the  next  section. 

The  time-domain  imaging  method  was  directly  imple¬ 
mented  using  Eq.  (4),  evaluated  using  straightforward  nu¬ 
merical  integration  over  all  incident- wave  and  measurement 
directions  employed.  The  synthetic  data  employed  was  sam¬ 
pled  at  rates  slightly  larger  than  the  Nyquist  frequency.  Be¬ 
fore  evaluation  of  the  argument  r  for  each  signal,  the  time- 
domain  waveforms  were  Fourier  interpolated  at  a  sampling 
rate  of  16  times  the  original  rate.  This  resampling,  as  well  as 
the  inverse  Hilbert  transform  from  Eq.  (4),  were  performed 
by  FFT.  Values  of  the  pressure  signals  at  the  time  r  were 
then  determined  using  linear  interpolation  between  samples 
of  the  resampled  waveforms. 

A  three-dimensional  point-spread  function  (PSF)  for  the 
present  time-domain  diffraction  tomography  method  is  il¬ 
lustrated  in  Fig.  1.  The  PSF  was  determined  by  recon¬ 
structing  an  ideal  point  scatterer  located  at  the  origin.  The 
time-domain  reconstruction  shows  a  dramatic  improvement 
over  the  single-frequency  reconstruction,  with  significant 
increases  in  both  the  point  resolution  (PSF  width  at  half¬ 
maximum  reduced  by  27%)  and  contrast  resolution  (first 
sidelobe  reduced  by  13  dB  and  second  sidelobe  reduced  by 
18  dB). 

Reconstructions  for  several  arbitrary  scattering  objects  are 
shown  in  Fig.  2.  All  of  these  reconstructions  were  performed 
using  synthetic  data  produced  by  the  fc-space  method  de¬ 
scribed  in  Ref.  [11].  Synthetic  scattering  data  were  com¬ 
puted  for  64  incident-wave  directions  and  256  measure- 
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Fig.  2.  Time-domain  reconstructions  from  full-wave  synthetic  data  for 
three  arbitrary  scattering  objects.  The  upper  row  shows  the  contrast  func¬ 
tion  7  for  each  object,  while  the  lower  row  shows  the  real  part  of  the 
reconstructed  contrast  7m*  Each  panel  shows  a  reconstruction  area  of 
5  mm  x  5  mm  using  a  linear  bipolar  gray  scale.  Left  to  right:  (a)  Cylinder, 
radius  2.5  mm,  with  an  internal  cylinder  of  radius  0.2  mm.  (b)  Cylinder, 
radius  2.5  mm,  with  random  internal  structure,  (c)  Tissue  structure,  with 
variable  sound  speed  and  density,  from  a  chest  wall  cross  section. 


ment  directions  in  each  case.  The  first  panel  shows  a  re¬ 
construction  of  a  cylinder  of  radius  2.5  mm  and  contrast 
7  =  -0.0295  with  an  internal  cylinder  of  radius  0.2  mm 
and  contrast  7  =  0.0632.  These  contrast  values  correspond, 
based  on  tissue  parameters  given  in  Ref.  [12],  to  the  sound- 
speed  contrasts  of  human  skeletal  muscle  for  the  outer  cylin¬ 
der  and  of  human  fat  for  the  inner  cylinder.  The  second  panel 
shows  a  reconstruction  of  a  2.5  mm-radius  cylinder  with  ran¬ 
dom  internal  structure.  The  third  reconstruction  shown  em¬ 
ployed  a  portion  of  a  chest  wall  tissue  map  from  Ref.  [13]. 
In  this  case,  the  synthetic  data  was  obtained  using  a  tissue 
model  that  incorporates  both  sound  speed  and  density  vari¬ 
ations,  so  that  the  actual  reconstructed  quantity  is  slightly 
different  from  7m  [1].  In  Fig.  2(c),  black  denotes  connec¬ 
tive  tissue,  dark  gray  denotes  muscle,  and  light  gray  denotes 
fat. 

The  real  part  of  each  reconstruction  in  Fig.  2  shows  good 
image  quality,  with  high  resolution  and  very  little  evidence 
of  artifacts.  Particularly  notable  is  the  accurately  detailed 
imaging  of  internal  structure  for  the  random  cylinder  and  the 
chest  wall  cross  section.  As  discussed  in  Ref.  [1],  the  den¬ 
sity  variations  present  in  the  chest  wall  cross  section  have 
not  greatly  affected  the  image  appearance;  there  is,  however, 
a  slight  edge  enhancement  at  boundaries  between  tissue  re¬ 
gions.  Also  notable  is  the  nearly-complete  absence  of  any 
artifacts  outside  the  scatterer  in  each  case;  this  result  indi¬ 
cates  that  high  contrast  resolution  has  been  achieved. 

Three-dimensional  reconstructions  of  a  homogeneous 
slab  with  sound  speed  contrast  7  =  0.01  and  dimensions 
1  mm  x  2  mm  x  3  mm,  are  shown  in  Fig.  3.  Synthetic  data 
was  computed  using  a  weak  scattering  approximation  for 
288  incident- wave  directions  and  1152  measurement  direc¬ 
tions,  each  evenly  spaced  in  the  angles  $  and  ©.  Isosurface 
renderings  of  the  real  part  of  the  reconstructed  7 m  are  shown 
for  the  surfaces  7 u  =  0.0025.  Consistent  with  the  point- 
spread  function  shown  in  Fig.  1,  the  time-domain  recon¬ 
struction  is  much  more  accurate  than  the  single-frequency 
reconstruction.  While  the  single-frequency  reconstruction 
shows  an  erroneously  rippled  surface,  the  time-domain  re¬ 
construction  is  smooth.  The  time-domain  reconstruction  is 
nearly  identical  to  the  original  object  except  for  some  round- 


Fig.  3.  Three-dimensional  reconstructions  of  a  uniform  slab  with  contrast 
7  =  0.01.  Each  reconstruction  shows  an  isosurface  rendering  of  the  surface 
7m  =  0.0025.  Left:  single-frequency  reconstruction.  Right:  time-domain 
reconstruction. 


Fig.  4.  Reconstructions  of  three  phantoms  from  measured  scattering  data. 
Each  panel  shows  an  area  of  9  mm  x  9  mm  using  a  bipolar  logarithmic 
scale  with  a  30  dB  dynamic  range.  Left  to  right:  (a)  Homogeneous  agar 
cylinder,  (b)  Agar  with  glass  spheres,  (c)  Agar  with  glass  spheres  and  three 
nylon  filaments. 

ing  of  the  sharp  edges  due  to  the  limited  high-frequency  con¬ 
tent  of  the  signal  employed.  The  length  scale  of  the  rounded 
edges  is  on  the  order  of  one-half  the  wavelength  of  the  high¬ 
est  frequency  in  the  pulse,  i.e.,  about  0.2  mm  for  the  —6  dB 
cutoff  of  3.25  MHz.  Notable  is  that  the  time-domain  method 
was  more  efficient  than  the  single-frequency  method  in  this 
case;  the  total  CPU  time  required  on  a  233  MHz  Pentium  n 
processor  was  100.0  CPU  min  for  the  time-domain  method 
and  233.4  CPU  min  for  the  single-frequency  method  (both 
computations  included  solution  of  the  applicable  linearized 
forward  problem  as  well  as  the  inverse  problem).  This  gain 
in  efficiency  was  possible  because  the  greatest  computa¬ 
tional  expense  occurred  in  the  “backpropagation”  of  the  sig¬ 
nals,  which  required  evaluation  of  complex  exponentials  for 
the  single-frequency  method,  but  only  linear  interpolation  of 
the  oversampled  farfield  pressure  waveforms  for  the  time- 
domain  method. 


IV.  Measurements 

The  practical  capability  of  the  time-domain  diffraction  to¬ 
mography  method  to  image  tissue-like  media  has  been  tested 
using  measured  scattering  data  for  three  tissue-mimicking 
phantoms,  each  of  diameter  6  mm.  Details  of  the  phan¬ 
tom  construction  and  measurement  procedure  are  given  in 
Ref.  [6]  and  briefly  summarized  here.  The  phantoms  are  pri¬ 
marily  composed  of  agar  (nominal  sound  speed  1510  m/s); 
one  is  homogeneous,  another  contains  tiny  (subresolution), 
randomly  distributed  glass  beads,  and  a  third  contains  three 
nylon  filaments  as  well  as  glass  beads.  Measurements  were 
made  using  a  ring  transducer  system  [14]  that  consists  of 
2048  elements,  each  of  which  can  be  used  independently 
as  a  transmitter  or  receiver.  This  fixed  transducer  config¬ 
uration  avoids  signal  degradation  from  phase  jitter  and  ex¬ 
cessive  scanning  time  associated  with  moving  transducers. 
The  control  electronics  associated  with  the  ring  transducer 
provide  the  capability  to  program  arbitrary  transmit  wave¬ 
forms.  The  element  pitch  is  0.23  mm,  less  than  one  half  of 
the  wavelength  at  the  nominal  center  frequency  of  2.5  MHz. 

Spatially-limited  plane  wave  pulses  were  transmitted 
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from  128  positions  equally  spaced  around  the  ring.  To  con¬ 
struct  the  spatially-limited  plane  waves,  a  4  mm-width  co¬ 
sine  rolloff  was  added  to  each  side  of  a  10  mm-width  uni¬ 
form  central  region  to  provide  a  smooth  transition  in  ampli¬ 
tude  and  reduce  wavefront  spreading.  A  backpropagation 
method  [15]  was  then  used  to  obtain  the  transmit  waveforms 
that  produced  the  desired  incident  wave. 

The  incident  field  (without  the  scattering  object)  and  the 
total  field  (with  the  scattering  object)  were  measured  around 
the  ring  for  each  incident  view.  To  compensate  for  sound 
speed  changes  due  to  water  temperature  variations,  the  back¬ 
ground  sound  speed  was  tracked  using  a  probe  beam  during 
the  measurement  of  both  the  incident  and  total  fields.  The 
sound  speed  in  the  background  was  estimated  from  knowl¬ 
edge  of  the  arrival  time  and  the  travel  distance  of  the  probe 
beam,  which  was  a  spatially  limited  plane  wave  directed 
to  the  side  of  the  phantom.  The  resulting  speed  estimate 
was  used  to  equalize  the  time  scale  of  all  waveforms.  A 
temperature-compensated  incident  field  p»(0 ,  a,  t)  was  sub¬ 
tracted  from  the  total  field  p(0 ,  a,  t)  to  obtain  the  scattered 
field  p8(a,0,t).  Finally,  wavefields  were  extrapolated  to 
128  measurement  positions  at  a  radius  of  7500  mm  by  an 
exact  spatio-temporal  transformation  [2] ,  [6] . 

Far-field  scattered  waveforms  for  each  incident-wave  di¬ 
rection  were  further  processed  by  a  deconvolution  opera¬ 
tion  [1]  that  compensated  for  transducer-dependent  varia¬ 
tions  in  the  incident  pulse.  The  result  for  each  incident- 
wave  direction  was  an  estimate  of  the  scattered  farfield  pres¬ 
sure  associated  with  an  ideal  incident  pulse  of  the  form 

f(t)  =  cos(27t/o£)  e~t2/(2<r2\  with  /0  =  2.25  MHz  and 
a  =  0.25  ps.  The  preprocessed  data  pfi(0,  a,  t)  were  then 
inverted  using  numerical  integration  of  Eq.  (4).  The  inver¬ 
sion  procedure  was  the  same  as  for  the  simulations  described 
above,  except  that  the  initial  sampling  rate  was  20  MHz  and 
that  signals  were  oversampled  to  80  MHz  by  Fourier  inter¬ 
polation. 

Reconstructions  for  the  three  phantoms  are  shown  in 
Fig.  4.  Each  panel  shows  good  reconstruction  quality  with 
a  uniform  background  and  high  point  and  contrast  resolu¬ 
tion  as  well  as  quantitative  accuracy  (similar  reconstruc¬ 
tions,  obtained  using  an  eigenfunction-based  inverse  scat¬ 
tering  method,  are  presented  in  Ref.  [6]).  The  subresolution 
glass  spheres  do  not  cause  speckle  as  in  pulse-echo  B-scan 
imaging,  but  instead  appear  as  slight  local  variations  in  con¬ 
trast  consistent  with  weak  point  scatterers.  Both  nylon  fila¬ 
ments  and  glass  spheres  appear  dark  because  higher  sound 
speed  corresponds  to  negative  contrast  7  as  defined  above. 
In  panel  (c),  reconstructions  of  the  nylon  wires  show  slight 
sidelobe  artifacts;  these  artifacts  could  be  removed  by  care¬ 
ful  choice  of  an  optimal  pulse  f(t)  in  the  preprocessing  of 
the  scattered  field  [6]. 

V.  Conclusions 

A  new  method  for  time-domain  ultrasound  diffraction  to¬ 
mography  has  been  presented  and  validated  using  synthetic 
and  measured  scattering  data.  The  method  provides  quanti¬ 
tative  images  of  sound  speed  variations  in  unknown  media. 
These  reconstructions  are  equivalent  to  multiple-frequency 
reconstructions  using  filtered  backpropagation,  but  can  be 
obtained  with  much  greater  efficiency.  The  time-domain  re¬ 
construction  algorithm  has  been  derived  as  a  simple  filtered 
delay-and-sum  operation,  closely  related  to  time-domain 
confocal  synthetic  aperture  imaging,  so  that  it  can  be  con¬ 
sidered  a  generalization  of  imaging  algorithms  employed  in 
current  clinical  instruments.  The  simplicity  of  the  imaging 


algorithm  allows  straightforward  addition  of  features  such 
as  time-gain  compensation  and  aberration  correction. 

Numerical  results  obtained  using  synthetic  and  measured 
data  show  that  the  time-domain  method  can  yield  signifi¬ 
cantly  higher  image  quality  (and,  in  some  cases,  also  greater 
efficiency)  than  single-frequency  diffraction  tomography. 
Quantitative  reconstructions,  obtained  using  signal  param¬ 
eters  comparable  to  those  for  present-day  clinical  instru¬ 
ments,  show  accurate  imaging  of  objects  with  simple  deter¬ 
ministic  structure,  random  internal  structure,  and  structure 
based  on  a  cross-sectional  tissue  model.  Reconstructions  of 
tissue-mimicking  phantoms  suggest  that  the  method  will  be 
useful  for  diagnostic  imaging  problems  such  as  the  detection 
and  characterization  of  lesions  in  ultrasonic  mammography. 
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A  finite-difference  time-domain  model  for  ultrasonic  pulse  propagation  through  soft  tissue  has  been 
extended  to  incorporate  absorption  effects  as  well  as  longitudinal-wave  propagation  in  cartilage  and 
bone.  This  extended  model  has  been  used  to  simulate  ultrasonic  propagation  through  anatomically 
detailed  representations  of  chest  wall  structure.  The  inhomogeneous  chest  wall  tissue  is  represented 
by  two-dimensional  maps  determined  by  staining  chest  wall  cross  sections  to  distinguish  between 
tissue  types,  digitally  scanning  the  stained  cross  sections,  and  mapping  each  pixel  of  the  scanned 
images  to  fat,  muscle,  connective  tissue,  cartilage,  or  bone.  Each  pixel  of  the  tissue  map  is  then 
assigned  a  sound  speed,  density,  and  absorption  value  determined  from  published  measurements  and 
assumed  to  be  representative  of  the  local  tissue  type.  Computational  results  for  energy  level 
fluctuations  and  arrival  time  fluctuations  show  qualitative  agreement  with  measurements  performed 
on  the  same  specimens,  but  show  significantly  less  waveform  distortion  than  measurements. 
Visualization  of  simulated  tissue-ultrasound  interactions  in  the  chest  wall  shows  possible 
mechanisms  for  image  aberration  in  echocardiography,  including  effects  associated  with  reflection 
and  diffraction  caused  by  rib  structures.  A  comparison  of  distortion  effects  for  varying  pulse  center 
frequencies  shows  that,  for  soft  tissue  paths  through  the  chest  wall,  energy  level  and  waveform 
distortion  increase  markedly  with  rising  ultrasonic  frequency  and  that  arrival-time  fluctuations 
increase  to  a  lesser  degree.  ©  1999  Acoustical  Society  of  America.  [S000 1-4966(99)032 12-9] 

PACS  numbers:  43.80.Qf,  43.80.Cs,  43.58.Ta,  43.20.Fn  [FD] 


INTRODUCTION 

Echocardiography  is  widely  employed  for  diagnosis  of 
cardiac  diseases  including  valvular  defects,  pericardial  effu¬ 
sion,  and  wall  motion  abnormalities.1-3  Commonly,  echocar¬ 
diography  is  performed  noninvasively  through  the  chest 
(transthoracic)  using  an  external  probe  placed  on  the  chest 
wall.  The  chest  wall,  however,  can  considerably  degrade  im¬ 
age  quality  because  acoustic  paths  between  the  skin  and 
heart  may  contain  ribs  and  cartilage  as  well  as  inhomoge¬ 
neous  muscle  and  fatty  tissue.  The  result  is  that  as  many  as 
10-30%  of  patients  cannot  be  successfully  imaged  with 
present  transthoracic  techniques.4  This  limitation  of  transtho¬ 
racic  echocardiography  has  led  to  the  development  of  transe¬ 
sophageal  echocardiography,  in  which  the  heart  is  imaged  by 
a  transducer  inserted  into  the  esophagus.1-4  Although  transe¬ 
sophageal  echocardiography  provides  superior  image  qual- 


a^Present  address:  Department  of  Meteorology,  The  Pennsylvania  State  Uni¬ 
versity,  University  Park,  PA  16802. 


ity,  resulting  in  high  diagnostic  sensitivity  and  specificity, 
the  invasiveness  of  the  procedure  is  accompanied  by  in¬ 
creased  risk.3-6  For  this  reason,  improvements  in  the  nonin- 
vasive  transthoracic  approach  are  desirable,  for  example,  by 
the  development  of  methods  to  compensate  for  image  degra¬ 
dation  caused  by  the  chest  wall. 

An  understanding  of  ultrasonic  aberration  produced  by 
the  chest  wall  is  important  to  the  development  of  appropriate 
compensation  methods  for  transthoracic  ultrasonic  imaging. 
Direct  measurements  of  ultrasonic  distortion  produced  by 
chest  wall  specimens7,8  have  been  helpful.  Results  reported 
in  Ref.  7  show  that  propagation  through  the  chest  wall  causes 
substantial  beam  distortion.  However,  that  study  did  not  dis¬ 
tinguish  the  effect  of  soft  tissue  from  effects  caused  by  rib 
structures.  In  Ref.  8,  a  detailed  study  of  distortion  caused  by 
soft  tissue  paths  indicates  that  soft  tissue  distortion  in  the 
chest  wall  is  substantially  less  than  the  corresponding  distor¬ 
tion  in  the  human  abdominal  wall.  However,  distortion 
caused  by  ribs  was  only  treated  qualitatively  in  the  latter 
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study  because  the  physical  mechanisms  of  rib-induced  dis¬ 
tortion  could  not  be  adequately  described  by  the  method  re¬ 
ported  there.  Although  a  model  of  ultrasound  propagation  in 
the  chest  wall  has  previously  been  described,9  that  model  is 
based  on  coarse  depictions  of  chest  wall  morphology  includ¬ 
ing  .  homogeneous  tissue  layers  and  evenly- spaced, 
uniformly- shaped  ribs.  These  previous  experiments  and 
simulations,  therefore,  have  left  gaps  in  the  current  knowl¬ 
edge  about  the  physical  causes  of  ultrasonic  wavefront  dis¬ 
tortion  caused  by  the  chest  wall. 

Recent  work  on  the  simulation  of  ultrasonic  pulse 
propagation10-12  has  provided  insight  about  the  wavefront 
distortion  caused  by  the  human  abdominal  wall.  Although 
these  studies  have  provided  specific  information  about  the 
relationships  between  soft  tissue  morphology  and  ultrasonic 
wavefront  distortion,  the  work  is  not  fully  applicable  to  dis¬ 
tortion  caused  by  the  human  chest  wall.  The  morphology  of 
chest  wall  soft  tissue  is  different  from  that  of  the  abdominal 
wall  in  ways  that  can  affect  ultrasonic  wavefront  distortion.8 
Furthermore,  imaging  through  the  chest  wall  is  complicated 
by  ribs  that  limit  the  usable  acoustic  window  size  and  cause 
scattering  and  reflection. 

The  study  reported  here  applies  quantitative  simulation 
methods,  similar  to  those  presented  in  Refs.  10  and  12,  to 
anatomically  detailed  chest  wall  models  that  include  the  ribs. 
Accurate  depiction  of  rib-ultrasound  interactions  requires 
not  only  representation  of  the  strong  reflections  associated 
with  sound  speed  and  density  contrast  between  ribs  and  soft 
tissue  (already  accurately  modeled  by  the  finite  difference 
method  of  Ref.  10),  but  also  modeling  of  the  strong  losses 
associated  with  propagation  through  bone  and  cartilage.  For 
this  reason,  the  finite-difference  method  described  in  Ref.  10 
has  been  extended  to  include  tissue-dependent  absorption. 
Quantitative  descriptions  of  the  distortion  caused  by  soft  tis¬ 
sues  are  obtained  by  statistical  analysis  of  simulated  distor¬ 
tion.  Visualizations  of  wavefronts  propagating  through  maps 
of  chest  cross  sections  provide  evidence  about  physical  rela¬ 
tionships  between  wavefront  distortion  and  the  morphology 
of  ribs  and  soft  tissue  structures  in  the  chest  wall.  Further 
insight  about  wavefront  distortion  mechanisms  is  provided 
by  a  comparison  of  distortion  results  for  incident  pulses  of 
different  center  frequencies. 


I-  THEORY 

Ultrasonic  pulse  propagation  through  the  human  chest 
wall  is  modeled  here  using  the  equations  of  motion  for  a 
fluid  of  variable  sound  speed,  density,  and  absorption.  The 
tissue  is  assumed  motionless  except  for  small  acoustic  per¬ 
turbations.  Absorption  is  included  using  an  adaptation  of  the 
Maxwell  solid  model,13  in  which  all  absorption  effects  are 
represented  by  a  single  relaxation  time.  This  assumption  re¬ 
sults  in  frequency-independent  absorption  characteristics. 
Equivalent  treatments  of  tissue-dependent  absorption  have 
been  employed  by  a  number  of  previous  models  for  ultra¬ 
sonic  propagation  in  biological  tissues.14-16  For  such  a  fluid, 
the  linearized  equations  of  mass  conservation,  momentum 
conservation,  and  state  can  be  combined  to  obtain  the  first- 
order,  two-dimensional,  coupled  propagation  equations, 


dp(x,y,t) 

Jt 


+  p(x,y)  c(x,y)2  V*v(x,y,r) 


=  -a(x,y)p(x,y,t), 
d\(x,y,t) 

p(x,y) - 37 - +  Vp(*,y  ,0  =  0. 


(1) 

(2) 


Here,  p(x,y,t)  is  the  acoustic  perturbation  in  fluid  pressure, 
v(x,y  ,0  is  the  vector  acoustic  particle  velocity,  p(x,y)  is  the 
ambient  density,  c(x,y)  is  the  ambient  sound  speed,  and 
a(x,y)  is  an  absorption  coefficient  that  is  equivalent  to  the 
inverse  of  a  spatially-dependent  relaxation  time  r(x,y). 

The  absorption  coefficient  a,  defined  as  a  real  quantity, 
is  related  to  the  energy  lost  per  unit  length  as  follows.  The 
propagation  equations  (1)  and  (2)  lead,  for  plane- wave 
propagation  of  the  form  p  =  el(kx~(dt )?  to  the  dispersion  rela¬ 
tion 


(3) 


where  k  is  the  complex  wavenumber,  co  is  the  (real)  radial 
frequency  2  irf,  and  c  is  the  (real)  sound  speed.  The  imagi¬ 
nary  part  of  the  wavenumber  k  is  the  absorption  in  nepers  per 
unit  length.  Thus,  the  absorption  parameter  a  can  be  ob¬ 
tained  by  a  numerical  solution  of  the  equation 


Im[fc]  = 


loss  (dB/length) 
20  log10(<?) 


(4) 


Solution  of  Eq.  (4)  results  in  wavenumbers  having  a  real  part 
that  differs  from  co/c .  Since  this  discrepancy  is  less  than  1% 
over  the  range  of  tissue  properties  employed  in  the  present 
study,  use  of  absorption  coefficients  computed  from  Eq.  (4) 
does  not  significantly  affect  propagation  characteristics  ex¬ 
cept  by  adding  the  specified  absorption. 

Equations  (1)  and  (2)  were  solved  numerically  using  the 
finite-difference  time-domain  (FDTD)  method  described  in 
Refs.  10  and  17.  This  method  is  a  two-step  MacCormack 
predictor- corrector  algorithm  that  is  fourth-order  accurate  in 
space  and  second-order  accurate  in  time.  The  computations 
employed  a  spatial  step  size  of  15  points  per  wavelength  at 
the  pulse  center  frequency  of  2.3  MHz.  Time  step  sizes  were 
computed  using  a  Courant-Friedrichs-Lewy  number  of 
0.25.  Further  details  on  this  class  of  finite  difference  algo¬ 
rithms  can  be  found  in  Refs.  18-20. 

The  initial  condition  was  chosen  to  model  the  experi¬ 
mental  configuration  in  Ref.  8,  in  which  a  spatially  broad, 
nearly  planar  wavefront  was  emitted  from  a  wideband, 
pulsed,  unfocused  source  far  from  the  tissue  layer.  The  initial 
wavefront  was  represented  in  the  present  simulation  as  a 
plane  wave  pulse  propagating  in  the  +y  direction: 


p(jc,y,0)  =  -sin[fc0(y-y0)]e  (y  ^o)2^), 


w(*,;y,0)  =  0, 


and 


v(x,y,0) 


Pp^y,  0) 

pc  ’ 


(5) 
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where  the  wavenumber  kQ  is  equal  to  Itt/q/c  for  a  center 
frequency  of  /0 ,  cr  is  the  Gaussian  parameter  of  the  pulse 
temporal  envelope,  and  u  and  v  are  the  x  and  y  components 
of  the  vector  acoustic  particle  velocity  v.  The  spatial  Gauss¬ 
ian  parameter  a  was  chosen  to  simulate  the  bandwidth  of  the 
pulse  used  in  the  experiments,  as  discussed  below  in  the 
Method  section. 

The  computational  configuration  is  analogous  to  that  de¬ 
scribed  in  Ref.  10.  The  domain  of  computation  is  two- 
dimensional,  with  the  y  direction  taken  to  be  parallel  to  the 
direction  of  propagation  and  the  x  direction  parallel  to  the 
initial  wavefront.  As  in  Ref.  10,  periodic  boundary  condi¬ 
tions  were  applied  on  the  domain  edges  that  were  parallel  to 
the  direction  of  propagation,  while  radiation  boundary  con¬ 
ditions  were  applied  on  the  edges  perpendicular  to  the  direc¬ 
tion  of  propagation. 

II.  METHOD 

This  study  employed  six  chest  wall  specimens  obtained 
during  the  autopsies  of  four  different  donors  between  79  and 
85  years  of  age  at  death.  One  specimen  (4L)  was  from  a 
white  female,  while  the  others  were  from  white  males.  After 
the  specimens  were  obtained,  they  were  stored  unfixed  at 
— 20  °C  and  thawed  when  needed  for  study.  Wavefront  dis¬ 
tortion  measurements  were  made  on  these  and  other  speci¬ 
mens  as  part  of  a  study  described  in  Ref.  8.  In  those  mea¬ 
surements,  2.3  MHz  ultrasonic  pulses  generated  by  a  0.5-in. 
piston  transducer  propagated  through  individual  chest  wall 
specimens  immersed  in  a  37  °C  water  bath  and  the  transmit¬ 
ted  pulses  were  measured  by  a  96-element  broadband  cardiac 
array  scanned  to  synthesize  a  two-dimensional  aperture.  Sta¬ 
tistics  describing  wavefront  distortion,  including  arrival  time 
fluctuations,  energy  level  fluctuations,  and  wave  shape  dis¬ 
tortion,  were  computed  for  the  measured  pulses. 

For  the  present  study,  six  of  the  previously  measured 
specimens  were  cut  into  ~7-mm  thick  cross  sections  using 
the  technique  described  in  Ref.  10.  The  slices  were  then 
fixed  and  stained  with  a  modified  Gomori’s  trichrome  stain 
according  to  the  procedure  detailed  in  Ref.  21,  so  that  tissue 
types  could  be  distinguished.  This  stain  colored  muscle  tis¬ 
sue  red  and  connective  tissue  blue  while  leaving  the  fat  its 
natural  color.  Calcified  tissue,  including  bone  and  cartilage 
in  the  current  specimens,  was  not  differentially  stained  by 
this  technique,  but  the  natural  contrast  between  bone,  carti¬ 
lage,  and  marrow  was  sufficient  to  allow  tissue  mapping. 
Full-color  300  d.p.i.  images  of  the  cross  sections  were  cre¬ 
ated  by  placing  each  stained  tissue  cross  section  directly  onto 
the  surface  of  a  flatbed  digital  scanner.  Image  editing  pack¬ 
ages  (Adobe  Photoshop,  Version  3.0,  and  the  Gnu  Image 
Manipulation  Program,  Version  1.0)  were  used  to  manually 
segment  the  cross  sectional  images,  i.e.,  to  map  the  images 
into  regions  that  corresponded  to  one  of  six  media.  The  me¬ 
dia  were  water  (representing  water  external  to  specimens  or 
blood  inside  blood  vessels),  fat  (including  subcutaneous  fat, 
fat  interlaced  within  muscle  layers,  and  marrow),  muscle, 
connective  tissue  (including  skin,  septa,  and  fasciae),  carti¬ 
lage,  and  bone  (including  cortical  bone  and  trabeculae  within 
cancellous  bone). 

3667  J.  Acoust.  Soc.  Am.,  Vol.  106,  No.  6,  December  1999 


The  nomenclature  employed  here  for  the  cross  sections 
corresponds  to  that  of  Ref.  8  for  the  whole  specimens  from 
which  the  cross  sections  were  taken;  each  cross  section  is 
identified  by  a  donor  number  together  with  “L”  or  “R”  to 
indicate  whether  the  corresponding  specimen  was  taken  from 
the  left  or  right  side  of  the  breastplate.  Additional  numbers 
were  used  in  Ref.  8  to  indicate  the  intercostal  space  used  in 
each  measurement;  here,  lower-case  letters  are  used  to  indi¬ 
cate  independent  acoustic  paths.  Wavefront  distortion  mea¬ 
surement  results  from  four  of  the  specimens  employed  here 
(4L,  5L,  7L,  and  7R)  were  reported  in  Ref.  8.  Distortion 
statistics  for  specimens  8L  and  8R  were  not  presented  in  Ref. 
8  because  of  limited  acoustic  windows.  No  new  measure¬ 
ments  were  made  for  the  present  study;  statistics  describing 
measured  distortion  are  taken  directly  from  Ref.  8. 

The  six  segmented  tissue  maps  are  shown  in  Fig.  1.  All 
of  the  cross  sections  contain  a  layer  of  septated  subcutaneous 
fat  below  the  skin.  Most  of  the  cross  sections  also  include  a 
layer  composed  primarily  of  the  major  pectoral  muscles  and 
their  connective  fasciae  above  the  ribs.  Between  the  ribs  are 
regions  of  muscle  (internal  intercostal  and  external  intercos¬ 
tal  groups)  interlaced  with  fat.  In  some  cases,  additional  thin 
layers  of  fat  between  muscle  layers  are  apparent.  Cross  sec¬ 
tions  4L  and  7R  are  cut  along  the  intercostal  spaces  parallel 
to  the  ribs,  so  that  in  each  a  wide  cross  section  of  soft  tissue 
appears.  Cross  sections  5L,  7L,  and  8L  are  cut  perpendicular 
to  the  ribs,  so  that  each  contains  soft-tissue  acoustic  paths 
with  width  equal  to  the  width  of  the  corresponding  intercos¬ 
tal  spaces.  Cross  section  8R  is  cut  perpendicular  to  the  ster¬ 
num  at  a  location  of  large  curvature  in  the  ribs,  so  that  the 
ribs  are  diagonally  sectioned.  Several  blood  vessels  appear  in 
cross  sections  4L,  7L,  7R,  and  8R;  the  largest  of  these  is  the 
internal  mammary  artery. 

The  basic  structure  of  the  cross  sections  is  consistent 
with  standard  descriptions  of  chest  wall  anatomy.22,23  Ribs 
appear  in  each  cross  section;  each  rib  is  composed  of  a  “cos¬ 
tal  cartilage”  near  the  sternum  (shown  in  most  of  the  cross 
sections  considered  here)  attached  to  a  “true  rib”  (composed 
primarily  of  cancellous  bone)  at  the  edge  farther  from  the 
sternum.  In  the  cross  sections  considered  here,  the  costal 
cartilages  are  primarily  composed  of  calcified  cartilage,  sur¬ 
rounded  by  a  thin  layer  of  cortical  bone  (solid,  dense  bone 
with  microscopic  porous  structure),  which  in  turn  is  sur¬ 
rounded  by  the  periosteum,  a  thin  membrane  of  connective 
tissue.  Cross  sections  7L  and  7R  also  appear  to  contain  a 
small  amount  of  cortical  bone  in  the  central  portion  of  the 
ribs.  This  phenomenon  may  be  associated  with  advanced  cal¬ 
cification  known  to  occur  in  aging  humans.24  Cancellous 
bone,  composed  of  thin  trabeculae  that  form  macroscopic 
cells  filled  with  marrow,  is  seen  in  all  the  ribs  of  cross  sec¬ 
tion  5L,  which  was  taken  at  a  distance  farther  from  the  ster¬ 
num  so  that  the  true  ribs,  rather  than  the  costal  cartilages, 
were  included  in  this  cross  section.  Some  cancellous  bone  is 
also  apparent  within  portions  of  the  ribs  of  cross  sections  4L 
and  8R.  In  each  case,  the  cancellous  bone  is  surrounded  by  a 
thin  layer  of  cortical  bone  and  by  the  periosteum.  A  portion 
of  the  sternum,  composed  of  cancellous  bone  surrounded  by 
cortical  bone,  is  visible  at  the  left  side  of  cross  section  4L. 

The  density  and  sound  speed  grids  needed  for  the  finite- 
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FIG.  1.  Chest  tissue  maps  used  in  simulations.  In  each  map,  blue  denotes  skin  and  connective  tissue,  cyan  denotes  fat,  puiple  denotes  muscle,  orange  denotes 
bone,  and  green  denotes  cartilage.  Blood  vessels  appear  as  small  water-filled  (white)  regions.  Simulated  apertures  are  indicated  using  lower-case  letters  for 
each  cross  section;  the  letters  correspond  to  the  acoustic  path  labels  used  throughout,  while  the  length  of  the  arrow  beneath  each  letter  corresponds  to  the  extent 
of  the  simulated  aperture.  Smaller  arrows  indicate  55-element  (11.60-mm)  apertures  while  large  arrows  indicate  68-element  (14.28-mm)  apertures. 


difference  computation  were  created  by  mapping  regions  of 
the  segmented  tissue  images  to  reference  density  and  sound 
speed  values  for  the  five  tissue  types  and  water.  The  water 
sound  speed  and  density  employed  are  those  of  pure  water  at 
body  temperature  (37.0  °C). 25,26  Sound  speeds  for  muscle 
and  fat  were  obtained  by  averaging  values  for  human  tissues 
given  in  Refs.  27  and  28.  A  representative  sound  speed  for 
connective  tissue  was  determined  using  an  empirical  formula 
relating  collagen  content  to  ultrasonic  sound  speed29  together 
with  a  measured  value  for  the  collagen  content  of  human 
skin.30  The  sound  speed  employed  for  bone  was  obtained 
from  an  average  of  values  reported  in  Ref.  31  for 
longitudinal-wave  propagation  in  human  cortical  bone.  The 
sound  speed  used  here  for  cartilage  is  that  given  in  Ref.  32  as 
quoted  in  Ref.  27.  Density  values  for  soft  tissues  were  deter¬ 
mined  from  Ref.  33  by  averaging  values  reported  for  adipose 
tissue,  skeletal  muscle,  and  skin,  respectively.  Density  values 
employed  for  bone  and  cartilage  are  average  values  from 
Ref.  31. 

Absorption  values  were  determined  from  attenuation 
measurements  summarized  in  Ref.  27  for  human  fat  at  37  °C, 
human  bicep  muscle  at  37  °C,  human  skin  at  40  °C,  human 
and  bovine  cartilage  at  23  °C,  and  human  skull  (temperature 
not  reported).  Attenuation  values  reported  at  other  ultrasonic 
frequencies  were  interpolated  (or,  for  the  skull  data,  extrapo¬ 
lated)  to  obtain  values  for  2.3  MHz  (corresponding  to  the 
pulse  center  frequency  employed  here  and  in  Ref.  8)  assum¬ 
ing  a  linear  dependence  of  attenuation  on  frequency.  This 
assumed  linear  dependence  is  a  simplifying  approximation; 
tissue  measurements  show  that  attenuation  varies  approxi¬ 
mately  as  where  j3  is  typically  between  0.9  and  1.5 


for  various  human  soft  tissues.34  The  absorption  for  water 
was  estimated  by  extrapolating  frequency-  and  temperature- 
dependent  absorption  values  summarized  in  Ref.  35  to  2.3 
MHz  and  37.0  °C.  The  values  of  tissue  parameters  employed 
in  the  present  study  are  given  in  Table  I. 

The  finite-difference  program  was  employed  to  compute 
propagation  of  a  plane  wave  pulse  through  each  scanned 
cross  section  from  the  skin  to  the  peritoneal  membrane, 
mimicking  the  propagation  path  employed  in  the  distortion 
measurements  of  Ref.  8.  The  spatial  step  size  of  the  finite- 
difference  grid  was  chosen  to  be  0.0442  mm,  or  1/15  wave¬ 
length  in  water  at  the  center  frequency  of  2.3  MHz.  The 
temporal  step  size  was  chosen  to  be  0.00725  fis,  for  an  op¬ 
timal  Courant-Friedrichs-Lewy  number  cAtfAx  of  0.25.20 
The  Gaussian  parameter  or  of  the  source  pulse  was  chosen  to 
be  0.4766  mm  in  accordance  with  the  experimentally  mea¬ 
sured  pulse  bandwidth  (for  pulses  transmitted  through  a  wa¬ 
ter  path)  of  1.2  MHz.  A  visual  comparison  confirmed  that  the 


TABLE  I.  Assumed  physical  properties  for  each  tissue  type  employed  in  the 
simulations. 


Tissue 

type 

Sound  speed 
(mm//is) 

Density 

(g/cc) 

Absorption 

(dB/mm) 

Water 

.  1.524 

0.993 

0.0007 

Fat 

1.478 

0.950 

0.12 

Muscle 

1.547 

1.050 

0.21 

Connective 

1.613 

1.120 

0.37 

Cartilage 

1.665 

1.098 

0.97 

Bone 

3.540 

1.990 

4.37 
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Simulated  pulse  closely  matched  the  measured  pulses  in 
shape  and  length. 

Each  simulation  was  performed  on  a  workstation  with 
128  MB  of  random-access  memory.  Finite-difference  grids 
on  the  order  of  1500X1000  points  were  employed.  At  each 
time  step,  the  wave  field  was  updated  on  a  grid  subset  chosen 
to  include  the  entire  support  of  the  acoustic  wave  but  to 
exclude  quiescent  regions.  The  entire  pressure  field  was 
saved  as  a  raster  image  at  intervals  of  0.725  /i s  for  later 
visualization.  The  computation  time  for  each  simulation  was 
on  the  order  of  five  hours.36 

Signals  were  recorded  for  8.62  (is  at  a  sampling  fre¬ 
quency  of  138  MHz  by  simulated  apertures  with  dimensions 
close  to  those  in  the  experimental  study  of  Ref.  8.  Positions 
of  all  simulated  apertures  employed  are  sketched  in  Fig.  1. 
The  simulation  of  receiving  elements  was  performed  by  in¬ 
tegrating  the  locally-computed  pressure  over  the  element 
pitch  of  0.21  mm.  For  cross  sections  cut  parallel  to  the  ribs, 
the  simulated  apertures  contained  68  elements  for  an  aperture 
width  of  14.28  mm.  For  cross  sections  cut  perpendicular  to 
the  ribs,  55  simulated  elements  were  used  to  form  11.55  mm 
apertures.  Element  directivity  effects  were  implicitly  incor¬ 
porated  by  the  integration  of  acoustic  fields  over  the  width  of 
each  element;  the  resulting  directivity  functions  correspond 
to  those  for  an  idealized  line  element  of  width  0.21  mm. 

A  one-dimensional  version  of  the  reference  waveform 
method10,37  was  used  to  calculate  the  arrival  time  of  the  pulse 
at  each  receiving  position  in  the  simulation  data.  In  this 
method,  the  relative  arrival  time  of  each  received  waveform 
is  computed  by  cross-correlation  with  a  reference  waveform. 
The  arrival  time  fluctuations  across  the  receiving  aperture  are 
then  calculated  by  subtracting  a  linear  fit  from  these  calcu¬ 
lated  arrival  times,  and  the  root-mean-square  value  of  these 
fluctuations  is  computed.  Energy  level  fluctuations  in  the 
data  were  calculated  by  summing  the  squared  amplitudes  of 
each  waveform  over  a  2.4- /is  window  that  isolated  the  main 
pulse,  converting  to  decibel  units,  and  subtracting  the  best 
linear  fit  from  the  resulting  values.  As  for  polynomial  fits 
previously  employed  in  wavefront  distortion  measurements,8 
the  purpose  of  the  linear  fit  removal  in  each  case  was  to 
compensate  for  gross  changes  in  tissue  thickness  across  the 
array.  Variations  in  pulse  shape  across  the  aperture  were 
evaluated  using  the  waveform  similarity  factor;37  this  quan¬ 
tity,  which  can  be  considered  a  generalized  cross-correlation 
coefficient,  has  a  maximum  of  unity  when  all  received  wave¬ 
forms  are  identically  shaped. 

To  test  the  frequency  dependence  of  chest  wall  wave- 
front  distortion,  propagation  through  eight  portions  of  speci¬ 
mens,  each  containing  only  soft  tissue,  was  also  computed 
for  wavefronts  having  center  frequencies  of  1.6  and  3.0 
MHz.  In  each  case,  the  initial  wavefront  was  chosen  to  have 
the  same  temporal  envelope  as  above.  The  absorption  coef¬ 
ficient  at  these  frequencies  for  each  tissue  type  was  extrapo¬ 
lated  from  the  value  employed  at  2.3  MHz  using  the  assump¬ 
tion  that  absorption  depended  linearly  on  the  center 
frequency.  The  spatial  and  temporal  sampling  rates  were  also 
varied  in  inverse  proportion  to  the  pulse  center  frequency. 
All  runs  were  otherwise  identical  in  configuration  and  pro¬ 
cessing  to  those  described  above. 
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TABLE  II.  Statistics  of  simulated  wavefront  distortion  caused  by  thirteen 
soft  tissue  paths  within  chest  wall  cross  sections.  The  “Path”  column  shows 
the  cross  section  label  and  aperture  letter  for  each  path;  these  labels  corre¬ 
spond  to  those  shown  in  Fig.  1.  The  statistics  shown  include  the  average 
specimen  thickness  for  the  tissue  path  considered,  rms  values  and  correla¬ 
tion  lengths  (CL)  of  the  arrival  time  fluctuations  (ATF)  and  the  energy  level 
fluctuations  (ELF),  the  waveform  similarity  factor  (WSF),  and  the  total  at¬ 
tenuation. 


Path 

Thickness 

(mm) 

ATF 

rms  CL 
(ns)  (mm) 

ELF 

rms  CL 
(dB)  (mm) 

Attenuation 
WSF  (dB) 

4L-c 

15.4 

32.0 

0.60 

1.98 

1.68 

0.981 

5.62 

4L-d 

12.7 

10.0 

2.58 

0.46 

1.23 

0.999 

4.08 

4L-e 

16.0 

10.0 

1.37 

1.61 

1.74 

0.998 

5.26 

4L-f 

17.0 

17.3 

2.48 

0.92 

1.61 

0.999 

5.33 

5L-a 

11.0 

11.6 

0.95 

1.51 

1.13 

0.991 

4.29 

5L-c 

15.0 

14.8 

1.03 

1.15 

1.19 

0.996 

5.01 

7L-a 

16.2 

16.8 

2.64 

0.95 

1.29 

0.999 

5.46 

7L-b 

14.9 

22.5 

2.66 

1.19 

1.61 

0.998 

4.91 

7R-c 

17.7 

17.4 

1.77 

2.52 

2.07 

0.997 

5.83 

7R-d 

21.0 

8.3 

1.10 

0.85 

1.79 

0.999 

7.07 

7R-e 

24.7 

13.7 

1.37 

1.06 

1.62 

0.997 

8.69 

8R-a 

23.8 

26.6 

1.78 

2.58 

1.40 

0.992 

7.76 

8R-b 

22.2 

29.9 

1.44 

1.95 

1.11 

0.989 

6.09 

Mean 

17.5 

17.8 

1.67 

1.44 

1.50 

0.995 

5.80 

St.  Dev. 

4.2 

7.8 

0.71 

0.66 

0.30 

0.005 

1.33 

III.  RESULTS 

Simulated  wavefront  distortion  results  for  13  soft  tissue 
paths  (i.e.,  paths  in  which  wavefront  distortion  was  not  sig¬ 
nificantly  influenced  by  the  ribs)  are  shown  in  Table  II. 
These  results  indicate  that  soft  tissue  paths  cause  a  wide 
range  of  wavefront  distortion  effects  depending  on  the  spe¬ 
cific  morphology  of  each  path.  For  instance,  path  7R-c 
causes  arrival  time  and  energy  level  fluctuations  that  are 
more  than  twice  the  magnitude  of  those  caused  by  the  adja¬ 
cent  path  7R-d.  This  difference  is  thought  to  arise  from  mor¬ 
phological  features,  including  muscle  tissue  with  interlaced 
fat  and  a  large  amount  of  connective  tissue,  of  the  tissue 
within  path  7R-c.  Also  notable  is  that  the  specimen  thickness 
does  not  closely  correspond  to  variations  in  distortion.  The 
largest  rms  arrival  time  fluctuation  and  lowest  waveform 
similarity  factor,  for  example,  are  caused  by  path  4L-c, 
which  has  an  average  thickness  less  than  the  mean  for  all  the 
tissue  paths. 

Wavefront  distortion  statistics  for  the  13  soft  tissue 
paths  are  graphically  summarized  in  Fig.  2  together  with 
corresponding  statistics  for  all  of  the  soft  tissue  measure¬ 
ments  reported  in  Ref.  8.  This  comparison  indicates  that 
wavefront  distortion  caused  by  soft  tissues  in  the  chest  wall 
simulations  is  comparable  to  measured  distortion.  Arrival 
time  fluctuations  and  energy  level  fluctuations  for  simulated 
distortion  are  slightly  less  than  measured  values,  but  mean 
values  of  both  fluctuations  for  the  simulations  fall  well 
within  one  standard  deviation  of  the  corresponding  mean 
fluctuation  for  the  measurements.  The  waveform  similarity 
factor,  however,  is  substantially  higher  for  simulations  than 
measurements,  indicating  that  simulated  waveforms  were 
distorted  considerably  less  than  measured  waveforms.  Corre¬ 
lation  lengths  for  the  simulated  distortions  are  somewhat  less 
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FIG.  2.  Summary  of  distortion  statistics  for  soft  tissue 
paths.  The  bar  chart  shows  mean  values  of  the  rms  ar¬ 
rival  time  fluctuations  (ATF),  rms  energy  level  fluctua¬ 
tions  (ELF),  correlation  lengths  (CL)  of  these  fluctua¬ 
tions,  and  waveform  similarity  factors  (WSF)  for  the 
simulations  performed  in  the  present  paper  and  the  ex¬ 
periments  reported  in  Ref.  8.  Error  bars  indicate  a  range 
of  plus  or  minus  one  standard  deviation  from  the  mean. 


than  measured  values.  However,  consistent  with  measure¬ 
ments,  the  mean  correlation  length  of  the  simulated  arrival 
time  fluctuations  is  greater  than  that  for  the  simulated  energy 
level  fluctuations. 

As  in  Ref.  8,  rib  structures  were  found  to  cause  much 
more  distortion  than  soft  tissue  alone.  The  varied  nature  of 
distortion  caused  by  rib  effects  is  illustrated  in  Fig.  3,  which 
shows  three  representative  sets  of  measured  signals  for  speci¬ 
men  8L.  These  measurements  were  made  during  the  study 
reported  in  Ref.  8.  The  first  panel  shows  96  adjacent  mea¬ 
sured  signals,  along  the  array  direction  (approximately  par¬ 
allel  to  the  ribs)  for  propagation  through  a  tissue  path  within 
an  intercostal  space.  The  signals  are  not  severely  distorted; 
secondary  arrivals  are  discernible,  but  are  of  lower  amplitude 
than  the  main  arrival.  The  second  panel  shows  96  measured 
signals  for  an  elevation  over  a  rib.  Here,  all  signals  are  se¬ 
verely  distorted.  Multiple  arrivals,  as  well  as  high-amplitude 
spatially-random  fluctuations,  are  seen.  The  third  panel 
shows  50  measured  signals  along  the  elevation  direction 
(perpendicular  to  the  ribs),  centered  over  the  soft  tissue  be¬ 
tween  the  ribs.  Here,  the  main  wavefront  is  curved  rather 
than  straight,  an  additional  arrival  behind  the  main  wavefront 
is  seen,  and  portions  of  the  signals  from  over  the  ribs  (at  both 
edges  of  the  panel)  are  advanced  relative  to  the  signals  from 
the  central  soft  tissue  region. 

The  present  simulations  allow  more  detailed  qualitative 
and  quantitative  investigation  of  rib  effects  than  were  pos¬ 
sible  from  the  previous  measurements.  Propagation  through 
two  rib-influenced  paths  is  illustrated  in  Figs.  4  and  5,  in 
which  computed  ultrasonic  pulses  are  superimposed  on  por¬ 
tions  of  the  tissue  maps  from  Fig.  1.  (Similar  visualizations 
of  propagation  through  soft  human  body  wall  tissue  were 
shown  in  Ref.  10.) 

Figure  4  shows  propagation  through  a  thin  rib,  com¬ 
posed  chiefly  of  cancellous  bone,  in  cross  section  5L  (corre¬ 
sponding  approximately  to  path  5L-b).  A  strong  reflection 
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FIG.  3.  Measured  waveforms  for  three  propagation  paths  in  specimen  8L. 
Each  panel  shows  received  waveforms  on  a  bipolar  logarithmic  gray  scale 
with  a  dynamic  range  of  40  dB.  The  horizontal  range  shown  in  each  panel  is 
20  mm  and  the  vertical  range  shown  is  6.4  fi s.  (a)  Tissue  path  between  two 
ribs,  in  azimuth  direction  (parallel  to  ribs),  (b)  Path  including  a  rib,  azimuth 
direction,  (c)  Tissue  path  including  intercostal  space  between  two  ribs,  el¬ 
evation  direction  (perpendicular  to  ribs). 
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FIG.  5.  Simulated  propagation  through  an  intercostal  space  in  cross  section  8L  (path  8L-b).  Panels  (a)-(d)  show  instantaneous  wavefields  at  successive 
intervals  of  3.62  jjls.  Each  panel  shows  an  area  that  spans  28.27  mm  horizontally  and  21.20  mm  vertically.  Wavefronts  are  shown  using  the  same  format  as 
in  Fig.  4. 


FIG.  4.  Simulated  propagation  through  the  central  rib  in  cross  section  5L  (path  5L-b).  Panels  (a)-(d)  show  instantaneous  acoustic  pressure 
intervals  of  2.17  jjls.  Each  panel  shows  an  area  that  spans  20.32  mm  horizontally  and  14.58  mm  vertically.  Logarithmically  compressed  wa 
on  a  bipolar  scale  with  black  representing  minimum  pressure,  white  representing  maximum  pressure,  and  a  dynamic  range  of  57  dB. 
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occurs  at  the  first  interface  between  bone  and  soft  tissue, 
removing  a  substantial  amount  of  energy  from  the  main 
wavefront.  The  small,  high-contrast  trabeculae  within  the  rib 
cause  considerable  scattering,  as  can  be  observed  in  panel  (b) 
of  Fig.  4.  The  scattering  causes  random  fluctuations  behind 
the  main  wavefront;  these  fluctuations  somewhat  resemble 
those  seen  in  the  measured  data  of  Fig.  3(b).  After  passing 
through  the  rib,  as  seen  in  panels  (c)  and  (d)  of  Fig.  4,  the 
central  portion  of  the  wavefront  shows  substantial  attenua¬ 
tion  and  distortion.  However,  the  average  arrival  time  of  the 
wavefront  is  not  greatly  changed  by  propagation  through  the 
rib,  but  is  advanced  by  only  about  one-half  period.  This  phe¬ 
nomenon  apparently  occurs  because  the  influence  of  the 
“slow”  marrow  (modeled  here  as  fat)  counteracts  the  influ¬ 
ence  of  the  “fast”  trabeculae.  Noteworthy  is  that  the  pre¬ 
dominant  ultrasonic  wavelength  has  increased  after  propaga¬ 
tion  through  the  rib,  so  that  the  effective  center  frequency  of 
the  wavefront  has  been  lowered.  Since  the  absorption  model 
used  in  the  present  study  includes  only  frequency- 
independent  absorption,  the  loss  of  short-wavelength  compo¬ 
nents  in  this  simulation  results  only  from  frequency- 
dependent  scattering  caused  by  the  trabeculae. 

Propagation  within  path  8L-b,  which  includes  two  larger 
ribs  and  the  corresponding  intercostal  space,  is  illustrated  in 
Fig.  5.  At  the  position  of  the  cross  section,  these  ribs  are 
composed  primarily  of  cartilage  and  surrounded  by  a  thin 
layer  of  cortical  bone.  Since  the  cartilage  and  bone  of  these 
ribs  are  modeled  as  homogeneous  structures,  small-scale 
scattering  within  these  tissues  did  not  occur  in  this  simula¬ 
tion.  Instead,  the  wavefront  is  reflected  from  interfaces  be¬ 
tween  cartilage,  bone,  and  soft  tissue. 

The  visualization  shown  in  Fig.  5  provides  physical  rea¬ 
sons  for  all  the  rib-related  distortion  phenomena  seen  in  the 
measured  data  of  Fig.  3(c).  The  wavefronts  propagating 
through  the  ribs  show  greater  attenuation  than  that  in  Fig.  4, 
both  because  of  the  high  absorption  of  the  ribs  and  because 
of  the  reflections  noted  above.  These  wavefronts  are  also 
advanced  relative  to  the  wavefront  propagating  through  the 
intercostal  space,  because  of  the  higher  sound  speed  of  both 
bone  and  cartilage.  The  wavefront  propagating  through  the 
intercostal  space  is  distorted  somewhat  by  the  inhomoge¬ 
neous  soft  tissue  path,  as  can  be  observed  in  panels  (b)  and 
(c).  However,  much  greater  distortion  results  from  interac¬ 
tion  between  the  wavefront  and  the  ribs.  A  rightward- 
propagating  reflection,  seen  in  panels  (b)  and  (c),  combines 
with  the  main  wavefront  in  panel  (d)  to  result  in  severe  dis¬ 
tortion  at  the  right  side  of  the  central  wavefront.  A  leftward- 
propagating  reflection  from  the  other  rib  is  also  apparent. 
Furthermore,  diffraction  from  the  edges  of  the  ribs  results  in 
large  curvature  of  the  soft  tissue  wavefront. 

Distortion  and  attenuation  statistics  for  a  variety  of 
simulations  employing  rib-influenced  paths  are  shown  in 
Table  HI.  Footnotes  in  Table  III  indicate  physical  causes  of 
distortion  present  within  each  path.  A  variety  of  distortion 
and  attenuation  mechanisms  are  illustrated.  Propagation 
through  small  intercostal  spaces  (paths  4L-a,  8L-b,  8L-f,  and 
7R-a)  causes  diffraction  effects  that  introduce  substantial 
curvature  into  the  wavefront,  as  seen  in  Fig.  5.  This  large- 
scale  wavefront  curvature  is  associated  with  large  arrival 
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TABLE  III.  Statistics  of  simulated  wavefront  distortion  caused  by  fourteen 
tissue  paths  including  rib  structures.  The  footnotes  associated  with  the  label 
for  each  path  indicate  morphological  features  and  physical  phenomena  that 
affected  the  wavefront  distortion  computed  for  that  path.  The  format  is 
analogous  to  that  in  Table  II. 


Path 

ATF 

Thickness  rms  CL 

(mm)  (ns)  (mm) 

ELF 

rms  CL 
(dB)  (mm) 

WSF 

Attenuation 

(dB) 

4L-aa,b-c'd 

21.0 

260.3 

3.00 

2.58 

2.72 

0.968 

15.33 

4L-bb'c 

17.6 

161.9 

1.90 

4.16 

1.49 

0.641 

43.35 

5L-bb 

14.2 

92.5 

0.69 

3.06 

1.92 

0.775 

26.87 

7L-cc,e 

17.8 

47.2 

1.58 

5.33 

2.04 

0.958 

19.66 

7R-aa,c,d 

30.4 

123.1 

2.12 

3.80 

1.78 

0.960 

16.57 

7R-bce 

24.3 

165.6 

2.71 

6.88 

2.07 

0.274 

43.06 

8L-ac 

25.3 

113.9 

1.18 

7.75 

2.29 

0.907 

32.44 

8L-ba,d 

22.8 

109.7 

2.05 

3.43 

1.22 

0.974 

10.28 

8L-cc 

28.8 

134.0 

2.75 

3.04 

1.57 

0.944 

40.47 

8L-dd 

23.6 

78.9 

0.64 

3.06 

1.55 

0.950 

6.78 

8L-ec 

26.4 

208.8 

1.91 

3.62 

1.50 

0.810 

44.27 

8L-fd 

28.5 

169.9 

1.79 

5.02 

1.95 

0.916 

10.70 

8L-gc 

27.6 

210.8 

1.40 

3.36 

1.35 

0.892 

44.22 

8R-cb,c 

24.9 

81.4 

2.08 

2.76 

1.25 

0.962 

44.32 

aSmall  intercostal  spaces. 
bCanceIIous  bone. 
cCortical  bone  and  cartilage. 
dStrong  rib  reflections. 
eCortical  bone  within  cartilage. 


time  fluctuation  values  although  the  wavefronts  generally  ap¬ 
pear  to  be  locally  smooth.  Interference  between  directly- 
transmitted  and  rib-reflected  wavefronts  (paths  4L-a,  8L-b, 
8L-d,  8L-f,  and  7R-a)  introduces  arrival  time,  energy  level, 
and  waveform  distortion  substantially  greater  than  that  for 
soft  tissue  paths  without  ribs.  Propagation  through  cancel¬ 
lous  bone  (paths  4L-a,  4L-b,  5L-b,  and  8R-c)  results  in  con¬ 
siderable  attenuation  and  large  waveform  distortion,  while 
propagation  through  cortical  bone  and  cartilage  (paths  4L-a, 
4L-b,  8L-a,  8L-c,  8L-e,  8L-g,  7L-c,  7R-a,  7R-b,  and  8R-c) 
results  in  even  larger  attenuation  but  smaller  distortion. 
Where  bone  is  embedded  within  cartilage  (paths  7L-c  and 
7R-b),  additional  scattering  also  occurs.  For  the  path  includ¬ 
ing  a  large  bone  inclusion  (path  7R-b),  this  scattering  results 
in  an  extremely  high  energy  level  and  waveform  distortion. 

Computed  frequency-dependent  wavefront  distortion 
statistics  are  summarized  in  Fig.  6.  Tissue  paths  used  for 
these  computations,  none  of  which  include  rib  structures,  are 
those  labeled  4L-d,  4L-f,  5L-a,  5L-c,  8R-a,  8R-b,  7L-a,  and 
7L-b  in  Fig.  1.  The  results  shown  in  Fig.  6  indicate  that 
arrival  time  fluctuations,  energy  level  fluctuations,  and  wave¬ 
form  distortion  all  become  more  severe  with  increasing  pulse 
frequency.  The  most  dramatic  change  is  in  the  energy  level 
distortion;  on  average,  the  rms  energy  level  fluctuations  for 
the  3.0-MHz  signals  are  2.3  times  those  for  the  1.6-MHz 
signals.  Correlation  lengths  of  both  arrival  time  and  energy 
level  fluctuations  decrease  with  frequency,  so  that  the  pre¬ 
dominant  length  scales  of  ultrasonic  wavefront  distortion  are 
seen  to  decrease  with  the  ultrasonic  wavelength.  As  with  the 
rms  distortion  statistics,  the  most  dramatic  frequency- 
dependent  change  is  in  the  energy  level  fluctuations.  Still, 
even  the  high-frequency  pulses  here  show  substantially 
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FIG.  6.  Summary  of  simulated  frequency-dependent 
distortion  results.  Mean  rms  arrival  time  fluctuations 
(ATF),  energy  level  fluctuations  (ELF),  correlation 
lengths  (CL)  of  these  fluctuations,  and  waveform  simi¬ 
larity  factors  (WSF)  are  shown  for  each  of  the  three 
pulse  frequencies  investigated.  Error  bars  indicate  a 
range  of  plus  or  minus  one  standard  deviation  from  the 
mean. 


smaller  distortion  than  that  previously  observed  in  experi- 

•  io~12  38 

ments  and  simulations  for  the  human  abdominal  wall. 

IV.  DISCUSSION 

As  with  earlier  simulations  of  propagation  through 
tissue,10,12  the  current  study  shows  qualitative  agreement 
with  measured  wavefront  distortion  results  for  similar 
specimens.8  However,  the  accuracy  of  the  present  model  is 
limited  by  simplifications  of  true  tissue  structure.  In  particu¬ 
lar,  the  computational  model  here  does  not  account  for  prop¬ 
erty  variations  within  tissue  types,  tissue  microstructure,  or 
three-dimensional  tissue  structure.  Each  of  these  simplifica¬ 
tions  limits  the  ability  of  the  present  model  to  precisely 
mimic  experimentally  measured  ultrasonic  wavefront  distor¬ 
tion.  These  limitations  are  discussed,  with  respect  to  soft 
tissues,  in  Ref.  10. 

The  modeling  of  ribs  adds  additional  complication.  In 
the  current  study,  individual  trabeculae  were  assumed  to  be 
composed  of  tissue  having  properties  identical  to  cortical 
bone,  an  assumption  known  as  Wolffs  hypothesis.39  The 
validity  of  this  hypothesis  has  been  questioned;40,41  however, 
measured  elastic  properties  of  individual  trabeculae  vary 
widely40’41  and  recent  work42  has  provided  support  for 
Wolffs  hypothesis.  Thus,  the  properties  employed  here  for 
trabecular  bone  can  be  regarded  as  reasonable  order-of- 
magnitude  estimates.  Likewise,  the  modeling  of  marrow  as 
fat  tissue  is  a  simplifying  assumption  that  may  have  limited 
validity,  although  available  data  suggest  that  the  density  and 
sound  speed  of  marrow  are  close  to  those  for  other  adipose 
tissues.31  In  addition,  the  present  model  for  cartilage  is  based 
on  measurements  of  normal  cartilage,  while  the  cartilage 
present  in  the  specimens  employed  here  was  calcified  due  to 
the  age  of  the  donors.  However,  density  measurements  made 
on  eight  representative  samples  of  calcified  cartilage  (two 
from  specimen  7R,  four  from  specimen  1R,8  and  two  from 
an  unused  specimen)  resulted  in  an  average  density  of 
0.00111  kg/m3,  which  is  different  by  only  1%  from  the  den¬ 


sity  assumed  here.  Since  sound  speed  in  calcified  tissue  has 

J  43  44 

been  empirically  shown  to  be  directly  related  to  density, 
this  small  change  in  density  suggests  that  the  acoustic  prop¬ 
erties  of  the  calcified  cartilage  in  our  specimens  is  close  to 
that  for  normal  cartilage. 

The  computations  reported  here  model  the  chest  wall  as 
a  fluid  of  variable  sound  speed,  density,  and  compressibility. 
This  model  implicitly  neglects  shear  wave  propagation.  The 
neglect  of  shear  waves  in  soft  tissues  is  believed  to  be  justi¬ 
fied  because  the  absorption  of  shear  waves  in  soft  tissues  is 
much  greater  than  absorption  of  longitudinal  waves  45,46  In 
calcified  tissues,  however,  significant  shear  waves  are  known 
to  be  generated  47,48  In  the  current  scattering  configuration, 
some  shear  waves  are  likely  generated  wherever  the  rib  sur¬ 
face  is  far  from  parallel  to  the  wavefront.  However,  since 
shear  wave  absorption  has  been  found  to  be  somewhat  larger 
than  longitudinal  wave  absorption  for  ultrasonic  propagation 
in  bone,47  the  significance  of  shear-wave  propagation  within 
bone  on  transmitted  ultrasonic  wavefronts  is  questionable. 
For  this  reason,  omission  of  nonlongitudinal  waves  in  the 
present  study,  as  in  another  computational  study  of  ultrasonic 
scattering  from  bone,49  is  believed  to  be  justified;  however, 
further  study  would  be  required  to  confirm  this  assumption. 

The  absence  of  frequency-dependent  absorption  is  a  pos¬ 
sible  source  of  error  in  the  present  estimates  of  total  tissue 
attenuation,  energy  level  fluctuations,  and  waveform  distor¬ 
tion.  However,  since  absorption  in  tissue  increases  approxi¬ 
mately  linearly  with  frequency,  lower  absorption  for  fre¬ 
quency  components  below  the  pulse  center  frequency  would 
nearly  cancel  higher  absorption  for  frequency  components 
above  the  center  frequency,  so  that  the  average  absorption 
incurred  by  a  wideband  pulse  should  still  be  computed  with 
fair  accuracy.  For  this  reason,  the  absence  of  frequency- 
dependent  absorption  in  the  calculations  reported  here  is  not 
considered  to  be  a  significant  source  of  error  in  the  computed 
attenuation  or  energy  level  fluctuation  curves.  Still,  the  in¬ 
clusion  of  frequency-dependent  absorption  would  result  in 
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additional  waveform  distortion  effects.  The  lack  of  this  effect 
is  a  likely  reason  for  the  lower  waveform  distortion  (higher 
waveform  similarity  factors)  obtained  from  simulations  as 
compared  to  measurements.  However,  the  absence  of 
frequency-dependent  absorption  effects  allowed  frequency- 
dependent  scattering  effects  to  be  clearly  quantified  sepa¬ 
rately  from  absorption  effects. 

Although  the  simulations  were  planned  to  match  the 
measurements  of  Ref.  8  closely,  a  number  of  differences  re¬ 
main.  The  most  important  of  these,  as  discussed  in  Ref.  10,  is 
that  the  simulations  were  performed  using  a  two-dimensional 
tissue  model  while  the  measurements  were  inherently  three- 
dimensional.  Other  differences  include  details  of  the  source 
waveform  and  wavefront  shape,  variations  in  the  specimen 
orientations  and  the  regions  interrogated,  and  variations  in 
the  distance  between  the  specimen  and  the  real  or  simulated 
receiving  aperture.  All  of  these  differences  could  contribute 
to  discrepancies  between  measurements  and  simulations. 

In  general,  most  of  the  simplifying  assumptions  in  the 
present  tissue  model  are  likely  to  result  in  underestimation  of 
wavefront  distortion  produced  by  the  human  chest  wall.  The 
incorporation  of  tissue  microstructure,  spatially-dependent 
acoustic  properties  for  each  tissue  type,  shear  wave  propaga¬ 
tion  in  bone  and  cartilage,  three-dimensional  propagation, 
and  frequency-dependent  absorption  could  all  result  in 
greater  spatial  and  temporal  variations  in  the  propagating 
acoustic  fields,  so  that  these  features  could  produce  simu¬ 
lated  distortion  with  characteristics  closer  to  measurements. 
For  this  reason,  distortion  statistics  computed  using  the 
present  tissue  model  should  be  interpreted  as  lower  limits  for 
the  statistics  of  distortion  occurring  in  real  chest  wall  tissue. 

Additionally,  some  of  the  discrepancy  between  simu¬ 
lated  and  measured  distortion  may  be  explained  by  the  non- 
uniform  characteristics  of  the  receiving  transducer  employed 
in  the  measurements.8  The  water-path  measurements  re¬ 
ported  in  Ref.  8  show  arrival  time  fluctuations  (mean  2.21 
ns)  and  energy  level  fluctuations  (mean  0.36  dB);  although 
small,  these  fluctuations  are  comparable  to  the  difference  be¬ 
tween  the  average  measured  and  simulated  fluctuations. 
Thus,  compensation  for  arrival  time  and  energy  level  fluc¬ 
tuations  due  to  transducer  irregularities  could  reduce  mea¬ 
sured  distortion  to  levels  closer  to  the  simulations.  Also,  the 
waveform  similarity  factor  for  water  path  measurements  was 
0.991, 8  which  indicates  greater  waveform  distortion  than  the 
average  value  of  0.995  computed  here  for  soft  tissue  paths. 
Thus,  compensation  of  the  measured  data  for  transducer 
impulse-response  variations  could  raise  the  measured  wave¬ 
form  similarity  factor  to  a  value  in  closer  agreement  with 
simulations. 

Previous  experimental  measurements  of  wavefront  dis¬ 
tortion  caused  by  the  human  chest  wall8  have  suggested  that 
distortion  caused  by  chest  wall  soft  tissues  is  less  severe  than 
that  caused  by  the  human  abdominal  wall.11’38  This  differ¬ 
ence  has  been  observed  to  occur  even  though  average  speci¬ 
men  thicknesses  were  comparable  in  chest  wall8  and  abdomi¬ 
nal  wall11,38  measurements.  The  present  results  provide 
support  for  these  results;  arrival  time  and  energy  level  dis¬ 
tortion  by  the  chest  wall  was  found  here  to  be  smaller  than 
that  produced  by  the  abdominal  wall  in  previous  simulation 


studies.10,12  For  the  simulations,  this  difference  may  be  paf- 
tially  explained  by  the  fact  that  the  chest  wall  specimens 
employed  here  are  thinner  on  average  (mean  thickness  17.5 
mm)  than  the  abdominal  wall  cross  sections  employed  in 
Refs.  10  and  12  (mean  thickness  26.7  mm).  Another  possible 
partial  explanation  is  that  the  pulse  center  frequency  em¬ 
ployed  in  abdominal  wall  measurements  and  simulations  was 
3.75  MHz,  significantly  higher  than  the  center  frequency  of 
2.3  MHz  for  the  chest  wall  measurements  and  simulations. 
Differences  in  pulse  frequency  and  specimen  thickness  may 
explain  the  discrepancy  in  energy  level  distortion  between 
the  abdominal  wall  and  chest  wall,  but  do  not  fully  explain 
the  discrepancy  in  arrival  time  distortion  results.  For  in¬ 
stance,  the  mean  arrival  time  and  energy  level  fluctuations 
per  unit  length  are  1.02  ns/mm  and  0.083  dB/mm  for  the 
present  study  vs  1.96  ns/mm  and  0.105  dB/mm  for  the  ab¬ 
dominal  wall  cross  sections  of  Ref.  10  and  12.  Arrival  time 
distortion  was  shown  here  to  increase  only  subtly  with  in¬ 
creasing  pulse  frequency,  so  that  this  discrepancy  in  arrival 
time  fluctuations  is  not  fully  explained  by  pulse  frequency 
differences.  However,  energy  level  fluctuations  increase 
markedly  with  frequency  for  chest  wall  tissue.  Thus,  for 
equal  ultrasonic  pulse  frequencies,  chest  wall  tissue  should 
cause  energy  level  distortion  per  unit  length  comparable  to 
that  caused  by  abdominal  wall  tissue. 

It  was  suggested  in  Ref.  8  that  chest  wall  morphology 
may  differ  from  abdominal  morphology  in  a  manner  that 
results  in  smaller  ultrasonic  wavefront  distortion.  The  cross 
sections  employed  here  can  be  compared  with  those  em¬ 
ployed  in  Refs.  10  and  12  to  evaluate  the  importance  of 
morphological  differences  between  chest  wall  and  abdominal 
wall  tissue.  One  difference  between  the  two  groups  of  cross 
sections  is  the  nature  of  the  subcutaneous  fat  layers.  The 
abdominal  wall  cross  sections  generally  contain  thicker  fat 
layers,  containing  many  more  lobular  structures  than  the 
chest  wall  cross  sections.  Since  the  high  contrast  between 
septa  and  fat  causes  substantial  ultrasonic  scattering,10-12  this 
morphological  difference  is  likely  to  result  in  lower  overall 
energy  level  and  waveform  distortion  for  chest  wall  tissue 
(although,  as  discussed  above,  the  energy  level  distortion  per 
unit  propagation  length  should  be  comparable).  Also,  the  ab¬ 
dominal  wall  and  chest  wall  cross  sections  have  a  markedly 
different  structure  within  the  muscle  layers  that  occur  below 
the  subcutaneous  fat.  The  abdominal  wall  cross  sections  have 
many  large-scale  features  due  to  aponeuroses  (interfaces  be¬ 
tween  muscle  groups,  composed  of  connective  tissue  and  fat) 
and  large  fatty  regions.  These  large-scale  features  cause  large 
wavefront  fluctuations  that  are  associated  with  large  rms  ar¬ 
rival  time  fluctuations.10,12  In  contrast,  muscle  layers  of  the 
chest  wall  cross  sections  considered  here  contain  primarily 
smaller-scale  structures  associated  with  small  islands  of  in¬ 
terlaced  fatty  tissue.  This  morphological  difference  may  re¬ 
sult  in  lower  large-scale  arrival  time  fluctuations  but  signifi¬ 
cant  energy  level  fluctuations  associated  with  scattering, 
consistent  with  the  differences  between  distortion  caused  by 
soft  tissues  in  the  abdominal  wall  and  the  chest  wall. 

The  present  results  for  the  frequency  dependence  of  dis¬ 
tortion  provide  further  insight  into  the  importance  of  scatter- 


3674  J.  Acoust.  Soc.  Am.,  Vol.  106,  No.  6,  December  1999 


Mast  et  a/.:  Simulation  of  pulse  propagation  3674 


irlg  effects  relative  to  large-scale  structure  in  wavefront  dis¬ 
tortion  caused  by  soft  tissues.  If  wavefront  distortion  in  the 
chest  wall  were  caused  only  by  large-scale  tissue  structures, 
the  distortion  would  be  roughly  independent  of  frequency, 
since  propagation  effects  are  independent  of  frequency  in  the 
geometric  acoustics  limit.  However,  distortion  caused  by 
scattering  effects  should  increase  with  the  pulse  frequency 
for  inhomogeneities  of  size  comparable  to  the  wavelength. 
Previous  simulation  and  experimental  studies10-12  on  distor¬ 
tion  caused  by  the  human  abdominal  wall  have  suggested 
that  energy  level  fluctuations  and  waveform  distortion  are 
generally  associated  with  scattering  effects,  while  arrival 
time  fluctuations  are  predominantly  caused  by  large-scale 
path  length  differences.  The  present  results,  while  consistent 
with  those  conclusions,  indicate  that  scattering  plays  a  role  in 
all  types  of  distortion  considered  here.  Since  energy  level 
fluctuations  and  waveform  similarity  factors  exhibit  more 
dramatic  increases  in  distortion  with  increasing  pulse  fre¬ 
quency,  the  present  results  suggest  that  scattering  is  of  pri¬ 
mary  importance  in  causing  energy  level  and  waveform  dis¬ 
tortion  and  of  secondary  importance  in  causing  arrival  time 
distortion. 

These  results  can  be  employed  to  evaluate  the  potential 
of  various  approaches  to  improve  echocardiographic  imag¬ 
ing.  Available  acoustic  windows  for  transthoracic  imaging 
are  severely  limited  by  the  presence  of  the  ribs,  so  that  image 
quality  cannot  be  significantly  improved  by  an  increase  of 
aperture  size.  The  present  results  also  indicate  that  use  of 
higher-frequency  probes  may  provide  less  benefit  than  ex¬ 
pected  because  of  frequency-dependent  scattering  in  the 
chest  wall. 

For  these  reasons,  aberration  correction  methods  are  po¬ 
tentially  important  in  transthoracic  echocardiography,  par¬ 
ticularly  for  higher-frequency  imaging.  The  frequency- 
dependent  distortion  results  reported  here  suggest  that 
distortion  models  employing  single  phase  screens  may  be  of 
some  benefit  for  aberration  correction  in  echocardiography 
through  soft  tissue  paths.  The  relatively  weak  dependence  of 
arrival  time  fluctuations  on  pulse  frequency  suggests  that  a 
large  portion  of  arrival  time  variations  are  caused  by  tissue 
structures  too  large  to  cause  significant  frequency-dependent 
scattering  effects.  Similar  conclusions  regarding  the  impor¬ 
tance  of  large-scale  structures  to  arrival  time  fluctuations 
have  also  been  drawn  from  results  presented  in  Refs.  10  and 
12. 

Still,  the  present  results,  like  those  from  earlier 
studies,10-12  suggest  that  single  phase  screens  will  not  pro¬ 
vide  complete  correction  for  distortion  caused  by  soft  tissues. 
In  particular,  methods  employing  single  phase  screens  will 
not  completely  remove  distortion  caused  by  scattering.  The 
sharp  increase  of  amplitude  and  waveform  distortion  with 
frequency,  as  well  as  the  moderate  increase  of  arrival  time 
distortion  with  frequency,  indicate  that  scattering  effects  be¬ 
come  much  more  important  to  ultrasonic  aberration  as  imag¬ 
ing  frequencies  increase.  Furthermore,  phase  screen  models 
do  not  inherently  account  for  distortion  caused  by  rib  struc¬ 
tures,  shown  here  to  produce  diffraction,  reflection,  and  scat¬ 
tering.  Thus,  any  attempted  correction  using  only  phase 


screen  models  is  likely  to  provide  little  improvement  in  the 
presence  of  strong  rib-induced  effects. 

Other  correction  models  that  incorporate  rib  structures 
may  provide  greater  image  improvements  for  the  distortions 
most  important  to  echocardiography.  Processing  wavefronts 
with  techniques  such  as  angular  spectrum  filtering  can  re¬ 
move  some  spurious  arrivals,50  although  such  computations 
may  be  difficult  to  incorporate  into  a  general  correction  al¬ 
gorithm.  Other  possible  methods  include  those  incorporating 
models  of  tissue  structure.  Models  incorporating  ray 
acoustics9  may  provide  improvement,  but  implicitly  neglect 
diffraction  and  scattering  effects,  so  that  aberration  correc¬ 
tion  would  be  incomplete,  particularly  for  small  intercostal 
spaces.  A  more  complete  aberration  correction  method  could 
employ  synthetic  focusing  using  full-wave  numerical  com¬ 
putation  of  acoustic  fields  within  sufficiently  accurate  models 
of  tissue  structure.  This  method  has  been  implemented, 
within  the  context  of  a  quantitative  frequency-domain  in¬ 
verse  scattering  method,  in  Ref.  51.  However,  the  results 
presented  here  indicate  that  distortion  caused  by  soft  tissue 
and  rib  structures  varies  widely  based  on  morphological 
variations  between  (and  within)  individuals.  Thus,  for  any 
general  correction  method  employing  models  of  tissue  struc¬ 
ture,  separate  models  of  tissue  structure  must  be  constructed 
for  each  region  of  interest. 

V.  CONCLUSIONS 

A  computational  study  of  ultrasonic  propagation  through 
the  chest  wall,  including  tissue-dependent  absorption  as  well 
as  detailed  anatomical  cross  sections,  has  been  presented.  For 
soft  tissue  paths,  computational  results  for  arrival  time  dis¬ 
tortion,  energy  level  distortion,  and  correlation  lengths  of 
these  distortions  are  comparable  to  those  reported  in  previ¬ 
ous  chest  wall  measurements.  Both  simulations  and  measure¬ 
ments  indicate  that  arrival  time  distortion  and  energy  level 
distortion  caused  by  soft  tissues  in  the  human  chest  wall  is 
smaller  than  that  caused  by  the  human  abdominal  wall.  Dif¬ 
ferences  in  morphology  between  the  abdominal  wall  and  the 
chest  wall  provide  a  probable  explanation  for  this  difference. 

Distortion  caused  by  rib  structures  is  much  more  severe 
than  that  caused  by  soft  tissues.  Reflections  and  diffraction 
from  rib  structures  complicate  wavefronts  that  travel  through 
soft  tissue  paths  adjacent  to  ribs  and  can  cause  arrival  time 
and  energy  level  fluctuations  much  greater  than  those  in¬ 
duced  by  soft  tissue  structures.  Wavefronts  propagating  di¬ 
rectly  through  rib  structures  are  attenuated  by  both  internal 
absorption  and  reflection  at  interfaces  between  bone,  carti¬ 
lage,  and  soft  tissue.  Internal  scattering  within  rib  structures 
causes  distortion  phenomena  that  include  severe  waveform 
and  energy  level  distortion,  additional  attenuation,  and  low¬ 
ering  of  the  effective  frequency  for  the  transmitted  pulse. 
The  strong  dependence  of  distortion  on  the  morphological 
details  of  rib  structures  presents  a  major  challenge  for  aber¬ 
ration  correction  in  echocardiography. 

Simulation  of  propagation  through  soft  tissue  paths  us¬ 
ing  three  different  pulse  frequencies  has  indicated  that  the 
distortion  types  investigated  here  have  different  frequency 
dependence.  Arrival  time  fluctuations  increase  subtly  with 
frequency,  while  energy  level  and  waveform  distortion  in- 
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crease  greatly.  Thus,  a  substantial  portion  of  arrival  time 
fluctuations  produced  by  the  chest  wall  may  be  explained  by 
large-scale  tissue  variations,  but  some  arrival  time  distortion 
and  most  energy  level  and  waveform  distortion  apparently 
result  from  scattering.  Thus,  correction  of  wavefront  distor¬ 
tion  caused  by  soft  tissues  should  become  both  more  impor¬ 
tant  and  more  challenging  as  pulse  frequencies  employed  in 
imaging  systems  are  increased. 
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Abstract  -  We  investigate  accuracy  of  existing  2D 
pseudospectral  and  k-space  formulations  for 
simulating  acoustic  propagation  in  tissue  or  model 
scattering  media.  They  are  intended  to  provide 
insight  into  tissue-ultrasound  interaction  and  a  “test 
bed”  for  aberration  correction  schemes  in  medical 
imaging.  Both  methods  employ  FFT’s  to  evaluate 
spatial  derivatives  to  high  accuracy  on  coarse  grids. 
The  primary  difference  lies  in  the  approach  to  time 
integration.  Scattering  in  large-scale,  2D, 
inhomogeneous  media  is  included.  We  compare 
simulations  against  analytical  solutions  to  illustrate 
spatial  and  temporal  discretization  required  for 
acceptable  solutions. 

INTRODUCTION 

The  medium  is  represented  by  a  uniform  Cartesian 
grid  where  pressure/stiffness  and  velocity/density 
are  unknowns/parameters  at  discrete  points. 
Spectral  operators  in  space  enable  accuracy  and 
computational  efficiency  in  very  large  models. 
However,  inhomogeneities  are  often  represented  as 
piecewise  constant  from  node  to  node,  rather  than 
smooth.  The  resulting  stairstep  can  produce 
spurious  diffractions  at  edges/comers,  inaccurate 
reflections  and  transmissions  at  interfaces  and  local 
Gibbs  phenomena,  by  approximating  derivatives  at 
a  material  discontinuity.  Thus,  the  efficiency 
permitted  by  coarse  spectral  grids  is  compromised 
by  the  need  to  resolve  interface  derivatives. 

For  example,  scattering  by  a  soft  cylinder  requires 
only  two  nodes  per  wavelength  inside  and  outside 
the  cylinder  for  accurate  propagation,  but 
significantly  more  nodes  per  wavelength  are 
necessary  to  reduce  interface  artifacts.  Interface 
artifacts  are  quantified  for  a  single  interface,  ID 


multilayer  models,  and  cylindrical  scatterers. 
Abdominal  wall  cross  sections  with  coarse  and  fine- 
scale  inhomogenities  are  used  to  explore  fidelity  of 
wave  propagation  versus  nodes  per  wavelength  and 
tissue  characteristic  lengths.  We  show  that  the 
existing  tools  are  useable  in  2D. 

The  pseudospectral  method  is  implemented  in  the 
SpectralFlex  code.  Kbench  implements  the  k-space 
method. 

PSEUDOSPECTRAL  AND  K-SPACE  METHODS 

The  pseudospectral  and  k-space  methods  were 
formulated  to  provide  efficient  high-accuracy 
solutions  to  long  range  wave  propagation  problems. 
In  fact,  they  debuted  during  the  same  year  [1,2].  We 
briefly  describe  the  two  methods  as  implemented  in 
[3,4],  highlighting  the  major  similarities  and 
differences. 

Both  use  FFT’s  to  evaluate  spatial  derivatives  to 
high  accuracy  on  coarse  grids.  The  primary 
difference  lies  in  their  respective  approaches  to  time 
integration.  Note  that  coarse  spatial  grids  provide 
the  primary  incentive  for  FFT  based  (or  any  high 
order)  method.  The  computational  burden  is  linear 
in  the  number  of  timesteps  per  cycle,  for  both  2D 
and  3D.  Including  the  timestep,  computational 
burden  is  proportional  to  the  number  of  Points  Per 
Wave  (PPW)3  in  2D  or  (PPW)4  in  3D. 

SpectralFlex  adopts  a  4*  order  staggered  Adams 
Bashforth  ABS4  time  integrator  [5].  Among 
general  purpose  integrators,  this  is  close  to  optimal 
for  the  current  applications  -  2-3  digits  of  accuracy 
for  a  wave  propagating  several  hundred  wavelengths 
on  the  coarsest  possible  grid.  The  stability  limit  for 
ABS4  in  2D  is  CFL  =  0.3.  The  CFL  number  is 
defined  as:  CFL  =  At/(Ax/c),  where  At  is  the 
timestep,  c  is  the  wavespeed  and  Ax  is  the  cell  size. 
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Accuracy  frequently  requires  a  smaller  timestep, 
say  CFL  =  0.1.  Kbench  implements  a  time 
integrator  in  k-space  based  on  the  exact  solution  for 
waves  propagating  in  a  homogeneous  medium  [4]. 
It  outperforms  the  general  purpose  ABS4  time 
integrator  for  weak  scatterers  in  a  homogeneous 
linear  acoustic  medium.  ABS4  becomes  more 
efficient  when  the  scattering  objects  have  a  larger 
impedance  contrast. 

LONG  RANGE  PROPAGATION 

To  illustrate  the  advantages  of  the  FFT  based 
approach  for  long  range  propagation,  we  propagate 
a  2.5  MHz  pulse  200  wavelengths  through  water 
using  both  SpectralFlex  and  PZFlex,  a  finite 
element  code  that  is  second  order  accurate  in  both 
space  and  time.  The  center  frequency  is  2.5  MHz, 
but  spectral  content  is  observable  up  to  5  MHz. 
Thus,  4  PPW  at  2.5  MHz  is  the  minimum  sampling 
capable  of  resolving  the  pulse. 

Figure  1  compares  exact,  PZFlex  and  SpectralFlex 
solutions.  SpectralFlex  used  4  PPW  at  CFL  =  0.1, 
whereas  PZFlex  used  20  PPW  at  CFL  =  0.8.  These 
discretizations  in  time  and  space  are  typical  of  those 
that  would  be  used  in  real  problems.  The 
SpectralFlex  signal  looks  good  and  can  be  made 
better  by  reducing  the  timestep.  The  PZFlex  signal 
is  delayed  in  time  and  badly  dispersed.  A  much 
finer  grid  is  required  to  achieve  reasonable 
accuracy.  Note  that  at  CFL  =1.,  PZFlex  becomes  a 
characteristic  method,  and  produces  exact  results, 
even  at  2  PPW.  Unfortunately,  this  only  works  for 
ID  linear  problems. 

Kbench  produces  exact  results  for  this  example 
because  the  time  integrator  is  based  on  the  exact 
solution  for  a  homogeneous  medium. 

DISCONTINUITIES 

Spectral  methods  compute  highly  accurate  spatial 
derivatives  of  smooth  fields.  Thus,  in  homogeneous 
regions,  2  cells  per  minimum  wavelength  (ie, 
highest  spatial  frequency)  suffice.  However,  at 
material  interfaces  both  the  pressure  and  velocity 
fields  should  exhibit  slope  discontinuities  as  given 
by  (1),  where  n  denotes  the  normal  direction  and  the 

dp+  p+  dp~ 

— —  (1) 

dn  p  dn 


superscript  defines  the  +  or  -  side  of  the  interface. 
The  velocity  field  likewise  exhibits  slope 
discontinuities  at  interfaces. 

Spectral  methods  enforce  smoothness, 
approximating  the  jumps  in  normal  derivatives  with 
steep  gradients  over  a  few  cells.  This 
approximation  is  quite  good  at  10-20  cells  per 
wavelength,  but  less  accurate  at  2  cells  per 
wavelength.  For  a  staggered  grid,  as  in 
SpectralFlex,  the  material  interfaces  coincide  with 
velocity  nodes,  so  we  average  the  density  at  these 
interface  points.  On  a  regular  grid,  all  the  nodes  lie 
away  from  interfaces,  so  no  averaging  is  necessary, 
but  the  accuracy  is  even  worse  than  for  the 
staggered  grid. 


Figure  1.  Long  range  pulse  propagation  through 
water. 


ID  versus  exact  solutions 

Table  1  summarizes  material  properties  used  for  the 
ID  benchmarks.  Figure  2  illustrates  the  reflection/ 
transmission  of  a  normally  incident  pulse  at  a 
water/fat  interface  as  modeled  by  SpectralFlex.  To 
plotting  accuracy,  the  transmitted  signals  appear 
exact  (because  it  has  much  larger  amplitude  than  the 
reflected  wave).  However,  the  error  in  the  reflected 
signal  is  readily  apparent  at  4  PPW,  and  barely 
visible  at  6  PPW. 

Figure  3  shows  results  for  a  water  /bone  interface. 
In  this  case,  errors  are  visible  in  both  the  reflected 
and  transmitted  signals  at  4  PPW.  In  both  codes,  the 
most  pathological  case  is  varying  density/constant 
stiffness.  Fortunately,  most  tissues  have  a  higher 
contrast  in  stiffness  than  density  [6],  so  this  worst 
case  is  seldom  encountered.  As  shown  in  Fig.  4 
(density=1000,  928  kg/m3)  errors  in  the  reflected 
wave  are  visible  even  at  12  PPW. 
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Table  1  -  Material  Properties  for  1-D  benchmarks 


Material 

Density  [kg/m3] 

Wavespeed  [m/sec] 

Water 

1000. 

1500. 

Fat 

928. 

1427. 

Conn 

1100. 

1537. 

Musi 

1041. 

1571.  i 

Livr 

1050. 

1577. 

Figure  2.  Reflected  pulse  at  a  water/fat  interface. 


Figure  3.  Reflected  pulse  at  a  water/bone  interface. 


Figure  4.  Reflected  pulse  from  worst  case  interface. 


The  next  benchmark  examines  propagation 
through  a  1-D  approximation  of  an  abdominal  cross 
section.  Material  parameters  are  again  given  in 
Table  1.  Slight  errors  in  the  transmitted  wave  are 
apparent  at  4  PPW  (Fig.  5),  but  not  at  8  PPW. 
Reflected  signals  (not  shown)  are  similar.  Figure  6 


illustrates  the  effect  of  coarse  non-conforming  grids. 
At  4.1  PPW,  cell  boundaries  are  misaligned  with 
actual  material  interfaces  by  up  to  Vi  cell.  This  is,  of 
course,  the  case  for  any  real  model  with 
discontinuous  material  properties.  Properties  are 
assigned  based  on  the  center  of  the  cell.  The  errors 
introduced  by  this  sampling  dwarf  all  others.  More 
will  be  said  about  this  in  a  later  section. 


Figure  5.  Pulse  transmitted  through  1-D  abdominal 
wall  model. 


Figure  6.  Pulse  transmitted  through  1-D 
approximation  of  abdominal  wall.  Non-conforming 
grid. 


Scattering  by  cylinders 

In  addition  to  the  numerical  errors  at  interfaces, 
approximations  are  introduced  by  the  stair-step 
representation  of  curved  surfaces.  To  quantify  these 
approximations,  we  consider  3  mm  radius  fat  and 
bone  cylinders  immersed  in  water  and  insonified  by 
the  usual  2.5  MHz  pulse.  We  compute  the  difference 
between  exact  and  numerical  signals  for  each 
timestep  at  128  locations  at  6  mm  radius,  and  equal 
spacing  in  theta.  We  use  the  L2  norm  of  this  matrix 
as  an  error  metric.  Figure  7a  shows  the  L2  error  vs 
PPW  for  kbench  and  SpectralFlex  at  CFL  =  0.2. 


1999  IEEE  Ultrasonics  Symposium 


(c)  1999  IEEE 


0-7803-5725-6/99/$1 0.00 


The  curves  are  similar,  though  kbench  is  slightly 
more  accurate.  For  the  larger  contrast  bone  case  in 
Figure  7b,  similar  trends  are  evident,  but  in  this  case 
SpectralFlex  is  more  accurate.  The  error  is  tending 
to  zero  as  the  PPW  increases.  The  rate  of 
convergence  is  not  quite  quadratic.  For  context,  Fig. 
12  shows  waveforms  for  L2  error  near  0.01. 


Table  2  -  Material  Properties  for  Cylinders 


Material 

Wavespeed  [m/sec] 

Density  [kg/m3] 

Water 

1524. 

993. 

Fat 

1478. 

950. 

Bone 

3540. 

1990. 

Figure  7c  illustrates  that  at  low  CFL,  the  error  due 
to  time  integration  tends  to  zero.  For  this  problem, 
kbench  permits  reasonable  accuracy  at  roughly 
double  the  SpectralFlex  timestep.  For  the  bone 
cylinder,  the  stability  limit  of  SpectralFlex  is  0.15 
(0.3  in  the  bone) ,  and  kbench  can  go  up  to  0.2. 


ULl 

a 


Cells  per  wavelength 
Fat  Cylinder,  CFL  =  0.2 


Cells  per  wavelength 
Bone  Cylinder,  CFL  =  0.1 


—  kbench 

—  Spectral 

-lex 

_ 

■ 

CFL  CFL 

Fat  Cylinder,  6  cells  per  wave  Bone  Cylinder,  6  cells  per  wave 

Figure  7.  Cylinder  benchmarks.  Convergence  with 
increasing  discretization. 


materials 


Tissue  examples 

Figure  8  shows  an  abdominal  wall  cross  section 
[7,8].  This  model  is  insonified  by  a  4.35  MHz  plane 
wave  pulse.  Figure  9  displays  typical  reflected  and 
transmitted  signals  computed  by  SpectralFlex  at  4,  8 
and  12  PPW.  The  grids  were  defined  such  that 
material  boundaries  always  lie  in  exactly  the  same 
place.  Again,  it  is  confirmed  that  even  the  coarse  4 
PPW  model  produces  fairly  accurate  results. 


Figure  9.  Transmitted  pulse  from  Abdominal  wall 
model. 


INTERFACE  TREATMENTS 

Given  that  the  largest  numerical  errors  in  the  FFT 
based  methods  stem  from  material  interfaces,  we 
look  at  several  interface  treatments  for  reducing 
those  errors. 

Jump  conditions 

One  possible  method  for  improving  the  accuracy  at 
interfaces  is  to  split  the  solution  into  smooth  and 
non-smooth  parts,  and  apply  the  spectral  method 
only  to  the  smooth  part.  The  idea  is  to  introduce 
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local  corrections  at  material  interfaces  that  enforce 
the  jump  conditions  exactly.  E.g.,  construct  low 
order  polynomials  over  the  cells  adjacent  to  the 
interface  that  have  zero  value  and  zero  slope  1  cell 
away,  and,  when  added  to  the  continuous  part 
satisfy  the  jump  condition  (1)  at  the  interface. 
Obviously,  the  correction  is  not  required  to  be  local, 
but  if  it  covers  more  than  1  cell,  the  algorithm  will 
become  much  more  complicated  for  multiple 
interfaces.  LaVeque  [9]  discusses  such  an  approach 
applied  to  finite  difference  models. 

Figure  10  compares  reflected  and  transmitted 
signals  for  coarse  models  of  an  interface  with  and 
without  the  jump  correction  for  interface  velocity. 
This  example  isolates  the  effects  of  density  changes 
in  that  only  the  density  is  discontinuous.  The  bulk 
modulus  is  continuous.  The  correction  term 
improves  the  computed  result,  but  not  to  the  level  of 
a  homogeneous  material.  A  similar  correction  could 
be  applied  to  the  discontinuity  in  the  velocity 
gradients.  However,  it  will  have  a  weaker  effect  on 
the  staggered  grid  since  the  leading  coefficients  are 
already  continuous. 


Figure  10.  Jump  treatment  applied  to  interfaces. 


Smoothing  (Bandlimitation') 

Another  approach  to  improving  accuracy  at 
discontinuities  is  to  smooth  or  bandlimit  the  model 
before  sampling.  This  removes  unresolvable  high 
spatial  frequencies  from  the  model  itself.  We  found 
that  perfect  bandlimitation  reduced  computed 
signals  too  much,  but  a  “halfband”  filter  improves 
accuracy.  The  halfband  filter  is  smooth  with  an 
amplitude  of  0.5  at  half  the  sampling  frequency. 
Figure  11  shows  direct  and  halfband  filtered  kbench 
models  of  a  3  mm  cylinder  using  the  same  number 
of  PPW.  The  corresponding  pressure  fields  are 
plotted  using  a  60  dB  bipolar  log  scale. 


a)  Unsmoothed  b)  Halfband  Filtered 

Figure  11.  Direct  sampled  &  bandlimited  cylinders. 
Models  and  pressure,  60  dB  bipolar  log  scale. 


The  staircase  representation  of  the  cylinder 
generates  diffracted  signals  at  each  comer  in  11a, 
but  these  have  disappeared  in  lib.  Figure  12 
shows  selected  waveforms  from  the  direct  and 
halfband  sampled  models.  The  late  time  diffractions 
have  been  removed,  and  overall  L2  error  was 
reduced  from  0.0155  to  0.0105.  This  exercise 
demonstrates  that  smoothing  can  be  useful. 
However,  there  are  some  practical  complications. 
The  current  procedure  computes  the  smoothed 
object  as  the  inverse  transform  of  the  object’s 
analytical  spectrum  multiplied  by  the  filter,  and  is 
thus  defined  only  for  objects  with  a  known 
analytical  spectrum.  The  extension  to  more  general 
models  defined  on  a  pixel  by  pixel  level  has  not  yet 
been  demonstrated.  Also,  continuous  variations  of 
material  properties  produce  a  large  number  of 
distinct  materials.  In  the  limit,  each  cell  of  the 
model  has  different  properties.  For  the  purely 
acoustic  case,  this  presents  little  difficulty,  but  when 
material  nonlinearity  or  viscoacoustic  damping  is 
added,  the  complexity  intensifies.  E.g.,  for  each 
wavespeed/damping  set,  an  optimization  problem 
must  be  solved  to  compute  the  appropriate 
relaxation  constants,  and  these  constants  must  be 
stored. 
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4  5  6  7 

time  (us) 


Figure  12.  Backscattered  signals  from  direct  (top) 
and  bandlimited  (bottom)  models. 

Note  that  this  procedure  adds  information  to  the 
model  compared  to  the  unsmoothed  case.  Because 
smoothing  is  applied  to  the  analytical  cylinder,  the 
continuous  variation  of  material  constants  provides 
a  richer  set  of  parameters  than  is  available  in  the 
unsmoothed  representation.  As  long  as  the  model  is 
known  to  higher  resolution  than  the  grid, 
information  will  be  added.  It  is  an  interesting 
question  whether  smoothing  would  be  beneficial  on 
a  grid  finer  than  the  pixel  by  pixel  model  definition. 
For  example,  the  UOR  tissue  cross  sections  [7,8]  are 
the  most  detailed  models  we  know  of.  These  are 
represented  as  piecewise  constant  with  a  pixel  size 
0.085  mm  (about  7  PPW  for  a  2.5  MHz  pulse).  For 
a  5  MHz  pulse,  the  coarsest  grid  would  have  finer 
resolution  than  the  model. 

Volume  averaging  of  material  constants  has  also 
been  shown  effective  [10].  This  adds  additional 
information  compared  to  the  unsmoothed  case,  and 
the  correction  is  more  local  than  smoothing. 
However,  the  practical  difficulties  are  the  same. 

As  a  last  resort,  increased  discretization  (brute 
force)  will  always  converge  to  an  accurate  solution. 
This  is  a  practical  solution  in  2D,  as  the  above  tissue 
examples  indicate. 


CONCLUSIONS 

Model  parameterization  is  a  critical  issue  and  puts 
all  the  above  results  in  practical  perspective.  As 
shown  above,  differences  in  material  constants  or 
interface  locations  cause  much  larger  differences  in 
reflected/transmitted  signals  than  any  numerical 
errors  in  the  FFT  based  methods.  For  gaining 
insight,  or  as  a  test-bed  for  aberration  correction 
schemes,  a  4  PPW  model  is  sufficient  at  frequencies 
of  2.5  MHz  or  greater.  Fine  grids  or  cell-by-cell 
representation  of  material  properties  are  needed  only 
for  more  accurate  rendition  of  model  geometry. 
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Abstract 


Large-scale  simulation  of  ultrasonic  pulse  propagation  in  inhomogeneous  tissue  is  important 
for  study  of  ultrasound-tissue  interaction  as  well  as  for  development  of  new  imaging  methods. 
Typical  scales  of  interest  span  hundreds  of  wavelengths;  most  current  two-dimensional  methods, 
such  as  finite-difference  and  finite-element  methods,  are  unable  to  compute  propagation  on  this 
scale  with  the  efficiency  needed  for  imaging  studies.  Furthermore,  for  most  available  methods  of 
simulating  ultrasonic  propagation,  large-scale  three-dimensional  computations  of  ultrasonic  scat¬ 
tering  are  infeasible.  Some  of  these  difficulties  have  been  overcome  by  previous  pseudospectral 
and  fc-space  methods,  which  allow  substantial  portions  of  the  necessary  computations  to  be  ex¬ 
ecuted  using  fast  Fourier  transforms.  This  paper  presents  a  simplified  derivation  of  the  fc-space 
method  for  a  medium  of  variable  sound  speed  and  density;  the  derivation  clearly  shows  the  rela¬ 
tionship  of  this  /c-space  method  to  both  past  fc-space  methods  and  pseudospectral  methods.  In  the 
present  method,  the  spatial  differential  equations  are  solved  by  a  simple  Fourier  transform  method 
and  temporal  iteration  is  performed  using  a  k-t  space  propagator.  The  temporal  iteration  proce¬ 
dure  is  shown  to  be  exact  for  homogeneous  media,  unconditionally  stable  for  “slow”  (c(x)  <  Co) 
media,  and  highly  accurate  for  general  weakly  scattering  media.  The  applicability  of  the  fc-space 
method  to  large-scale  soft-tissue  modeling  is  shown  by  simulating  two-dimensional  propagation 
of  an  incident  plane  wave  through  several  tissue-mimicking  cylinders  as  well  as  a  model  chest  wall 
cross  section.  A  three-dimensional  implementation  of  the  A; -space  method  is  also  employed  for  the 
example  problem  of  propagation  through  a  tissue-mimicking  sphere.  Numerical  results  indicate 
that  the  fc-space  method  is  accurate  for  large-scale  soft-tissue  computations,  with  much  greater 
efficiency  than  that  of  an  analogous  leapfrog  pseudospectral  method  or  a  2-4  finite  difference  time- 
domain  method.  However,  numerical  results  also  indicate  that  the  /c- space  method  is  less  accurate 
than  the  finite-difference  method  for  a  high-contrast  scatterer  with  bone-like  properties,  although 
qualitative  results  can  still  be  obtained  by  the  /c-space  method  with  high  efficiency.  Possible  exten¬ 
sions  to  the  method,  including  representation  of  absorption  effects,  absorbing  boundary  conditions, 
elastic-wave  propagation,  and  acoustic  nonlinearity,  are  discussed. 
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I.  Introduction 


Computation  of  a  scattered  acoustic  field,  given  an  incident  wavefield  and  complete  specification 
of  an  inhomogeneous  medium,  is  known  as  the  forward  scattering  problem.  Numerical  solution  of 
the  forward  scattering  problem  is  central  to  many  aspects  of  ultrasonic  imaging,  including  inverse 
scattering  methods,  numerical  studies  of  wavefront  distortion,  and  development  of  new  methods 
for  adaptive  focusing.  Most  methods  for  numerical  solution  of  the  forward  scattering  problem 
fall  into  one  of  three  categories:  finite-difference  methods,  finite-element  methods,  and  spectral 
methods. 

Finite-difference  and  finite-element  methods  are  known  as  local  because  the  wave  propagation 
equations  of  interest  are  solved  at  each  point  based  only  on  conditions  at  nearby  points.  In  contrast, 
spectral  methods  such  as  the  fc-space  method  [1]— [7]  and  the  pseudospectral  approach  [8]— [14] 
are  called  global  because  information  from  the  entire  wavefield  is  employed  to  solve  the  wave 
propagation  equations  at  each  point.  In  part  because  of  their  global  nature,  spectral  methods  can 
be  more  accurate  than  local  methods — for  instance,  pseudospectral  methods  applied  to  periodic 
problems  have  been  shown  to  be  equivalent  to  finite-difference  methods  of  infinite  order  [12]. 

Spectral  methods  also  have  considerable  advantages  for  large-scale  forward  solvers  because 
the  required  storage  and  the  number  of  operations  per  iteration  can  be  dramatically  reduced  com¬ 
pared  to  local  methods.  This  advantage  occurs  principally  because  spectral  methods  can  allow 
computations  to  be  performed  on  coarser  grids  while  maintaining  accuracy.  For  example,  finite- 
element  methods  and  high-order  finite-difference  methods  typically  require  grid  spacings  on  the 
order  of  ten  points  per  minimum  wavelength,  while  second-order  finite-difference  methods  can 
require  twenty  points  per  wavelength  [10].  Spectral  methods,  in  theory,  require  only  two  points 
per  wavelength  (spatial  Nyquist  sampling),  although  for  computations  of  propagation  in  inhomo¬ 
geneous  media,  greater  accuracy  is  achieved  with  grid  spacings  on  the  order  of  four  points  per 
wavelength  [10, 11, 14]. 

This  report  addresses  the  problem  of  large-scale  ultrasonic  wave  propagation  in  biological  me¬ 
dia  such  as  human  tissue.  For  problems  of  interest  in  medical  ultrasound,  domain  sizes  can  often 
exceed  the  capabilities  of  conventional  forward  solvers.  For  example,  one  computation  of  realis¬ 
tic  scale  would  be  the  simulated  propagation  of  a  pulse  with  an  upper  bandwidth  limit  of  5  MHz 
in  a  volume  of  dimensions  30  mm  on  each  side  and  a  nominal  sound  speed  of  1.5  mm//js,  so 
that  the  minimum  wavelength  is  0.3  mm.  For  this  computation,  a  second-order  finite-difference 
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method  (using  twenty  points  per  wavelength)  would  require  a  three-dimensional  grid  containing 
8x  109  nodes,  a  finite-element  or  fourth-order  finite-difference  method  (using  ten  points  per  wave¬ 
length)  would  require  1  x  109  nodes,  and  a  spectral  method  (using  four  points  per  wavelength) 
would  require  6.4  x  107  nodes.  Since  a  grid  of  6.4  x  107  single-precision  complex  numbers  re¬ 
quires  storage  of  512  megabytes,  only  spectral  methods  are  feasible  for  realistic  three-dimensional 
computations  on  present-day  computers  that  typically  have  a  maximum  random-access  memory 
storage  of  several  gigabytes.  The  efficiency  provided  by  fast  Fourier  transform  implementations  of 
spectral  algorithms  is  a  further  reason  why  spectral  methods  are  a  practical  approach  to  large-scale 
and  three-dimensional  computations  of  ultrasonic  wave  propagation. 

Previous  spectral  approaches  have  included  pseudospectral  methods,  in  which  spatial  deriva¬ 
tives  are  evaluated  globally  by  Fourier  transformation  and  wavefields  are  advanced  in  time  using 
various  numerical  integration  techniques  [8]— [14].  This  method  has  provided  high  accuracy  in 
many  cases;  however,  temporal  iteration  techniques  that  provide  good  accuracy  for  large-scale 
models  typically  require  small  time  steps,  significant  additional  computations,  or  storage  of  wave- 
fields  from  additional  time  steps  [13],  so  that  the  efficiency  advantages  of  the  pseudospectral  ap¬ 
proach  are  less  than  might  first  be  expected.  The  k- space  family  of  methods  [1]— [7]  can  overcome 
this  problem  by  providing  explicit  temporal  propagators  related  to  the  Green’s  function  for  wave 
propagation  in  k-t  (spatial  frequency  and  time)  space. 

The  present  paper  presents  a  simplified  derivation  of  the  A;-space  method  using  a  differential 
representation  of  the  wave  propagation  equations.  The  spatial  part  of  the  wave  propagation  equa¬ 
tions  is  solved  by  Fourier  transformation  in  a  manner  analogous  to  past  pseudospectral  methods; 
this  derivation  is  shown  to  be  theoretically  equivalent  to  previous  integral  formulations  of  the  k- 
space  method.  Temporal  iteration  is  performed  using  a  k-t  space  propagator  [2],  which  is  shown 
to  be  exact  for  homogeneous  media,  and  in  general  to  provide  much  greater  accuracy  and  stability 
than  leapfrog  iteration  (in  which  temporal  derivatives  are  evaluated  using  second-order-accurate 
finite  differences)  without  significant  additional  computation  or  storage  requirements.  Thus,  the 
fc-space  method  provides  spatial  and  temporal  accuracy  ideal  for  large-scale  models  of  acoustic 
propagation  in  weakly  scattering  media. 

Below,  a  derivation  of  the  fc-space  method  is  presented  for  propagation  in  fluid  media  with 
spatially-dependent  sound  speed  and  density.  For  several  canonical  forward  problems  relevant 
to  ultrasonic  imaging,  the  accuracy  and  efficiency  of  the  /c-space  method  is  compared  to  a  pseu¬ 
dospectral  method  employing  leapfrog  iteration  and  to  a  2-4  finite  difference  time-domain  method. 
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The  A;-space  and  finite-difference  methods  are  also  used  in  an  example  computation  for  a  large- 
scale  two-dimensional  tissue  model.  Another  example  computation  illustrates  the  efficiency  of  the 
fc-space  method  for  three-dimensional  scattering  computations.  Possible  extensions  of  the  present 
fc-space  method,  including  multiple  relaxation  effects  for  absorption,  absorbing  boundary  condi¬ 
tions,  inclusion  of  elastic  and  nonlinear  acoustic  effects,  and  parallelization,  are  discussed. 
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II.  Theory 

A.  Derivation  of  the  A;-space  method 


The  /c-space  method  for  solving  the  acoustic  scattering  problem  is  briefly  derived  below.  The 
derivation  is  simpler  than  those  previously  published,  and  also  provides  some  new  insight  regarding 
the  remarkable  accuracy  and  stability  characteristics  of  the  /c-space  method. 

The  method  is  applicable  to  large-scale  modeling  of  linear  ultrasonic  propagation  in  soft  tis¬ 
sues,  which  are  modeled  here  as  fluid  media  with  spatially-dependent  sound  speed  and  density.  Al¬ 
though  the  /c-space  method  described  below  can  be  extended  to  include  absorption  effects,  acoustic 
nonlinearity,  and  shear-wave  propagation,  these  effects  are  neglected  in  this  derivation  for  simplic¬ 
ity. 

For  a  fluid  medium  with  spatially-dependent  sound  speed  and  density,  the  linear  acoustic  wave 
equation  is  [15] 

V- 

where  p(x,  t)  is  the  acoustic  perturbation  in  pressure,  p(x)  is  the  spatially-dependent  density,  and 
c(x)  is  the  spatially-dependent  sound  speed. 

By  defining  the  normalized  wavefield  /(x,  t)  =  p(x,  t) /^/p(x),  as  performed  in  a  number  of 
previous  studies  (e.g.,  Refs.  [16]  and  [17]),  the  first-order  derivative  term  can  be  eliminated  from 
Eq.  (1).  Details  of  this  step  are  given  in  Ref.  [6].  When  the  wavefield  is  also  split  into  incident  and 
scattered  parts,  such  that  /(x,  t)  =  /j(x,  t)  +  /s(x,  t),  a  wave  equation  for  the  scattered  field  can 


1  _  ,  A  1  d2p(x,t) 

—  Vp(x,t)J  -p(x)c(x)2  ft2  ~  ’ 


then  be  written 


V7s(x,*)--o 


1  d2fs{x,t ) 


x,t)  + 


d2v{x,ty 


v  jsk~,vj  $  at2  c?0vv  ’ '  1  dt2  r 

The  terms  on  the  right-hand  side  are  effective  sources  associated  with  density  and  sound  speed 


variations,  defined  as 

q(x,  t )  =  cl  y/p(x)  V2  ^1/ yfpip^j  f  (x,  t) 


(3) 


and 

V(X’*)  =  (c^  -1)/(x><)-  (4) 

The  incident  wavefield  fi(x,  t )  is  required  to  satisfy  the  usual  wave  equation  without  any  source 
terms  (i.e.,  the  D’Alembertian  operator  from  the  left-hand  side  of  Eq.  (2),  applied  to  fi(x,t), 
is  equal  to  zero).  Thus,  the  total  wavefield  f(x,  t )  also  satisfies  Eq.  (2)  identically,  so  that  the 
numerical  algorithm  developed  below  for  the  scattered  field  is  equally  applicable  to  the  total  field. 
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With  the  additional  definition  of  an  auxiliaiy  field  tu(x,  t)  =  fa(x,  t )  +  u(x,  t),  Eq.  (2)  can  be 
rewritten  in  k- space  as  the  coupled  set  of  equations 


=  M)2  [V(M  -  ^M)]  -  Q(M), 


m2 

V(M)  =  F 


[(- 


c(x)s 

Co 


Q( k,  t)  =  Co  F 


vfev2 


[/i(x,f)+w(x,f)] 

(vp5o)  ^ + ^  -  ^x>  ^ 


(5) 

(6) 

(7) 


where  F  denotes  spatial  Fourier  transformation  and  capital  letters  indicate  spatially  Fourier  trans¬ 
formed  quantities. 

For  each  point  in  fc-space,  Eq.  (5)  represents  an  independent  ordinary  differential  equation 
equivalent  to  the  standard  simple-harmonic  oscillator  equation  with  the  source  terms  (cok)2  V  and 
-Q.  This  ordinary  differential  equation  can  be  discretized  in  several  ways.  For  instance,  a  second- 
order  accurate  finite-difference  representation  of  the  second-order  time  derivative  allows  Eq.  (5)  to 
be  written  as 

Q(k,t) 


W(k,  t  +  At)  -  2  W( k,  t)  +  W{ k,  t-At)m  (cokAtf 


V(k,t)-W(k,t)~ 


(co  ky 


(8) 


where  Af  is  the  time  step.  This  is  known  as  “leapfrog”  iteration;  use  of  Eq.  (8)  in  the  cur¬ 
rent  method  is  analogous  to  commonly  used  pseudospectral  approaches  [13,  14].  (Although  in¬ 
creased  accuracy  can  be  achieved  by  higher-order  methods  such  as  fourth-order  Adams-Bashforth 
or  Adams-Moulton  iteration,  these  methods  have  the  disadvantage  of  requiring  storage  of  the  entire 
computational  grid  for  additional  time  steps  [12, 13].) 

A  more  accurate  form  of  the  temporal  iterator  is  obtained  using  a  nonstandard  finite  difference 
approach.  For  the  homogeneous  simple  harmonic  oscillator  equation,  an  exact  discretization  is 
known  [18].  (That  is,  for  any  temporal  and  spatial  step  sizes,  the  discrete  difference  equations  yield 
exactly  the  same  solutions  as  the  continuous  differential  equations.  A  similar  exact  discretization 
for  the  linear  part  of  the  Korteweg-de  Vries  equation  was  presented  in  Ref.  [19].)  Use  of  this 
nonstandard  discretization  leads  to  the  following  discrete  form  of  Eq.  (5): 


W{k,t+At)-2W{k,t)+W(k,t-At)  =  4sin2 


V(k,t)-W(k,t) 


ofcty 

(cok)2  _ 


•  (9) 


Because  the  discretization  employed  is  exact  for  the  simple  harmonic  oscillator  equation, 
Eq.  (9)  is  exactly  equivalent  to  the  differential  equation  (5)  for  the  case  of  a  homogeneous  medium 
[i.e.,  V(k,t)  =  Q(k,  t )  =  0].  Numerical  results  shown  below  indicate  that  high  accuracy  is  also 
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achieved  for  weakly  scattering  media,  in  which  case  V (k,  t)  <€.W (k,  t)  and  Q(k,  t)  <  W (k,  t). 
The  present  discretization  method  is  equivalent  to  that  employed  by  Bojarski  (the  form  given  in 
Ref.  [2]  follows  after  some  trigonometric  manipulation);  however,  previous  derivations  of  this 
method  have  been  based  on  approximations  to  an  integral  representation  of  Eq.  (5)  [2,  6].  It  may 
also  be  noted  that  Eqs.  (8)  and  (9)  are  equivalent  in  the  limit  of  small  At.  However,  results  shown 
below  indicate  that,  for  weakly  scattering  media,  use  of  the  k-t  propagator  (9)  provides  much 
greater  accuracy  for  larger  time  steps. 

In  numerical  implementation  of  the  fc-space  algorithm,  Eq.  (5)  is  used  to  advance  the  auxiliary 
field  W(k,  t)  in  time.  Equations  (6)  and  (7)  represent  updates  of  the  effective  scattering  sources 
v  and  q  and  their  spatial  Fourier  transformation  to  yield  the  A;- space  effective  sources  V  and  Q. 
Notable  is  that  the  effective  source  v  is  directly  proportional  to  the  square  of  the  sound  speed 
variation  of  the  medium,  while  the  effective  source  q  is  directly  proportional  to  the  Laplacian  of 
1  /\fp{x).  Thus,  for  a  piecewise-constant  inhomogeneous  medium,  v  may  be  non-zero  everywhere 
while  q  is  nonzero  (and  singular)  only  on  borders  between  regions. 

The  present  k- space  algorithm  can  now  be  summarized  as  follows: 

1 .  Set  any  initial  conditions  for  w  (x,  t )  and  spatially  Fourier  transform  (by  FFT)  to  obtain  initial 
conditions  for  W (k,  t). 

2.  Define  the  incident  wave  /j(x,  t)  on  the  entire  grid  (/j(x,  t )  can  be  identically  zero). 

3.  Compute  u(x,  t)  and  transform  to  obtain  V (k,  t)  [Eq.  (6)]. 

4.  Compute  q(x,  t)  and  transform  to  obtain  Q( k,  t)  [Eq.  (7)]. 

5.  Evaluate  W(k,t  +  At)  [Eq.  (9)]  and  inverse  transform  to  obtain  w(x,  f  +  At). 

6.  Set  t  ->•  t  +  At  and  go  to  step  (2). 

This  method  requires  three  fast  Fourier  transform  operations  per  time  step  (one  each  for  steps  3, 4, 
and  5  of  the  algorithm  enumerated  above). 

Also  notable  is  that  the  algorithm  is  directly  applicable  to  one-dimensional,  two-dimensional, 
and  three-dimensional  propagation.  This  is  possible  because  the  k-t  space  Green’s  function  has  an 
identical  form  for  any  number  of  spatial  dimensions  [2].  For  example,  to  implement  the  present 
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methods  for  two-dimensional  computations,  the  algorithm  outlined  above  is  simply  employed  us¬ 
ing  two-dimensional  Fourier  transforms.  The  three-dimensional  version  of  the  algorithm  is  for¬ 
mally  identical,  but  with  three-dimensional  Fourier  transforms. 

To  distinguish  between  the  standard  leapfrog  iteration  method  and  the  improved  method  used 
here,  the  following  nomenclature  is  used  in  the  present  paper.  The  above  algorithm  employing 
Eq.  (9)  is  referred  to  as  a  k- space  method,  while  the  corresponding  algorithm  employing  Eq.  (8) 
for  temporal  iteration  is  referred  to  as  a  leapfrog  pseudospectral  method.  This  nomenclature  is 
used  because  the  algorithm  employing  Eq.  (9)  is  mathematically  equivalent  to  an  extended  form  of 
Bojarski’s  A;-space  method  [2]  cast  in  terms  of  differential  equations  rather  than  integral  equations. 
The  algorithm  employing  Eq.  (8)  is  referred  to  as  pseudospectral  because  it  is  mathematically 
equivalent  to  a  conventional  “method  of  lines”  pseudospectral  algorithm  with  leapfrog  iteration 
[12].  (A  conventional  pseudospectral  algorithm  of  this  form  would  employ  the  spatial  inverse 
Fourier  transform  of  Eq.  (8)  for  temporal  iteration.) 

B.  Temporal  and  spatial  sampling  criteria 

To  determine  the  usable  range  of  spatial  and  temporal  sampling  rates  for  the  present  k- space 
method,  limits  on  the  stability  and  accuracy  of  the  method  can  be  evaluated  analytically. 

The  stability  of  the  fc-space  and  leapfrog  pseudospectral  methods  derived  above  can  be  evalu¬ 
ated  using  standard,  linear  von  Neumann  stability  analysis  [20].  In  this  technique,  the  difference 
equations  that  comprise  Eqs.  (8)  and  (9)  are  applied  to  a  test  function 

Wtest  (k,  nAt)  =  d  (k)n  ip  (k) ,  (10) 

where  ip( k)  is  a  spatial-frequency  domain  eigenmode  and  #(k)  is  a  temporal  amplification  factor. 
If  a  difference  equation  admits  solutions  with  |$(k)|  >  1  for  any  vector  wavenumber  k,  errors  may 
grow  exponentially  with  time  and  the  solution  is  thus  unstable.  If  |i9(k)|  <  1  for  all  wavenum¬ 
bers,  then  the  solution  is  numerically  stable.  For  simplicity,  the  present  stability  computation  is 
performed  in  the  absence  of  density  variations;  the  incident  wave  /j(x,  t)  is  assumed  (without  loss 
of  generality)  to  be  zero.  To  obtain  limiting  stability  criteria,  the  worst-case  sound  speed  inhomo¬ 
geneity  c(x)  =  Cmax  is  also  assumed. 

Application  of  this  technique  to  Eq.  (8),  which  represents  a  leapfrog  pseudospectral  approach, 
yields  a  quadratic  equation  for  #(k).  The  resulting  stability  condition  is 

Cmax^max At  ^  2,  (11) 
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where  (w  is  the  maximum  sound  speed  in  the  region  of  computation,  kmax  =  ir/Ax  is  the  max¬ 
imum  wavenumber  in  the  discrete  Fourier  transforms  used  to  compute  VF(k,t),  and  At  and  Ax, 
respectively,  are  the  temporal  and  spatial  steps  employed.  Using  the  standard  definition  for  a 
Courant-Friedrichs-Lewy  (CFL)  number  [21],  the  stability  condition 


CFL  =  ^ 

A.X  7T  Cmax 


(12) 


is  obtained  for  the  leapfrog  pseudospectral  method  represented  by  Eq.  (8). 

Application  of  the  same  analysis  to  the  k- space  iterator  of  Eq.  (9)  yields  the  stability  condition 


.  ttCFL 
sin  — - —  < 
2 


Co 

Qnax 


(13) 


This  condition  has  the  remarkable  result  that,  for  media  with  c(x)  <  Co  everywhere,  the  linear 
numerical  stability  of  the  k- space  method  is  unconditional.  However,  for  any  medium,  an  upper 
limit  on  the  time  step  still  arises  from  the  requirement  of  sampling  at  the  Nyquist  rate:  that  is,  the 
time  step  should  be  sufficiently  small  to  allow  two  samples  per  period  for  the  highest-frequency 
component  of  the  computed  field.  Thus,  the  temporal  sampling  criterion  can  be  written 


At  < 


7T 


2/ir 


Cmax^n 


Ax 

Cmax 


(14) 


or  simply  CFL  <  co/cmax.  The  stability  criterion  (13)  is  met  whenever  the  Nyquist  sampling 
criterion  (14)  is  met;  thus,  the  Nyquist  sampling  criterion  is  more  restrictive. 

For  the  spatial  discretization,  a  Nyquist  criterion  based  on  the  maximum  spatial  frequency 
kmax  =  k  I  Ax  is  met  for  any  step  size  Ax.  However,  the  inhomogeneous  medium  will  be  inaccu¬ 
rately  represented  (aliased)  if  its  Fourier  transform  has  significant  spatial-frequency  components 
beyond  kmax.  Aliasing  is  a  particular  problem  when  the  medium  contains  discontinuties  (which 
correspond  to  infinite  spatial-frequency  content);  removal  of  errors  associated  with  discontinuities 
is  discussed  in  the  following  section. 


C.  Effects  of  Discontinuities 

The  Fourier  transforms  performed  in  the  present  k- space  algorithm  can  lead  to  numerical  artifacts 
(related  to  the  Gibbs  phenomenon)  when  the  inhomogeneous  medium  contains  discontinuities  in 
sound  speed  or  density.  To  avoid  such  artifacts,  the  scattering  object  can  be  spatially  filtered  to 
smooth  any  discontinuities.  That  is,  the  spatially-dependent  sound  speed  c(x)  and  density  p(x) 
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can  be  replaced  by  filtered  functions  of  the  form 


«filtered(x)=F-1[£/(k)<^(k)], 


(15) 


in  which  the  Fourier  transform  U{ k)  of  the  function  it(x)  is  multiplied  by  a  low-pass  spatial- 
frequency  filter  </>(k).  The  function  U(k)  should  be  represented  as  accurately  as  possible;  for 
example,  exact  Fourier  transforms  of  simply-shaped  inhomogeneities  can  be  used  when  available. 
Below,  the  exact  Fourier  transform  of  a  two-dimensional  disk  is  employed  for  filtered  representa¬ 
tions  of  an  infinite  cylinder. 

In  the  present  study,  the  filter  employed  is  the  half-band  filter  [22] 


= 


'  1, 

k  f(k/kmax  -  1/2), 


k/kmax  <  1/2, 

1/2  <  k/kmax  <  3/2, 


(16) 


where 

f(0)  =  \  +  ^  cos(tt0)  -  ^  cos(3t r0)  (17) 

and  k  is  the  magnitude  of  the  spatial-frequency  vector  k. 

This  filter  defines  a  smoothly  tapered  window  that  causes  no  attenuation  of  spatial  frequencies 
below  kmax/2  and  drops  to  half  amplitude  (-6  dB)  at  the  spatial  frequency  kmax.  Zero  ampli¬ 
tude  is  reached  at  the  spatial  frequency  3/2  kmax,  which  exceeds  the  spatial-frequency  range  of  the 
discrete  Fourier  transforms  employed  in  the  fc-space  algorithm,  so  that  aliasing  error  is  not  elim¬ 
inated  by  the  half-band  filter.  However,  a  strict  bandlimiting  filter  was  found  to  cause  excessive 
loss  of  high-spatial-frequency  components  in  the  medium,  so  that  scattering  amplitude  near  the 
backscatter  direction  was  greatly  reduced.  The  half-band  filter  of  Eq.  (16)  was  found  to  greatly 
reduce  Gibbs  phenomenon  artifacts  while  maintaining  enough  high-spatial-frequency  components 
of  inhomogeneities  to  provide  accurate  backscatter  results. 

For  inhomogeneous  media,  exact  Fourier  transforms  are  not  generally  available.  However, 
artifacts  associated  with  discontinuities  can  still  be  removed  by  the  methods  given  above.  For 
example,  a  finely  sampled  representation  of  the  medium  could  be  filtered  using  Eq.  (15)  and  then 
decimated  to  the  desired  spatial  step  size. 
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III.  Numerical  Methods 


Numerical  implementation  of  the  k- space  algorithm  was  accomplished  using  the  algorithm  de¬ 
scribed  above.  The  normalized  incident  wave  /j(x,  t)  was  defined  as  a  plane  wave  with  Gaussian 
temporal  shape, 

fi(x,  t )  =  p(x)"2  sm(w0r) e-r2/(2<r2),  (18) 

where  r  is  the  retarded  time  r  =  t  -  (x  -  x0)/co  and  x0  is  the  initial  central  position  of  the 
wavefield.  This  incident  wave  was  implicitly  specified  using  initial  conditions  (as  for  the  incident 
plane  wave  in  Ref.  [23])  rather  than  explicitly  updated  at  each  time  step.  Boundary  conditions 
were  implicitly  periodic  at  each  edge  of  the  computational  domain,  due  to  the  inherent  periodicity 
of  the  fast  Fourier  transforms  employed. 

Wavefields  were  computed  on  two-dimensional  grids  large  enough  to  avoid  influence  of 
“wraparound”  error  within  the  temporal  window  of  interest.  All  k- space  computations  were  per¬ 
formed  on  square  grids  of  size  N  by  N.  Prior  to  execution  of  the  main  computation  loop,  the  Lapla- 
cian  occurring  in  Eq.  (7)  was  evaluated  using  second-order  accurate,  centered  finite-difference 
representations  of  the  second  derivative  in  each  direction.  Within  the  main  computational  loop, 
all  spatial  derivatives  were  evaluated  by  Fourier  transformation,  implemented  using  a  fast  Fourier 
transform  (FFT)  algorithm  [24].  For  maximum  FFT  efficiency,  grid  sizes  N  were  chosen  to  be 
integers  with  prime  factors  no  higher  than  3. 

To  reduce  any  spatial  anistropy  associated  with  the  rectangular  grid  shape,  the  spatial-frequency 
time-domain  wavefield  W{k,  t  +  At)  was  windowed  using  the  radially  symmetric  window 

(j>{ k)  =  H(kmax  -  k)  (19) 

before  inversion  to  yield  w(x,  t  + At).  (That  is,  within  step  5  in  the  algorithm  enumerated  above.) 
In  Eq.  (19),  H  is,  as  before,  the  Heaviside  step  function,  kmsx  is  the  maximum  wavenumber  mag¬ 
nitude  (equal  to  n/Ax,  since  the  spatial-frequency  range  sampled  extends  from  -tt/Ax  to  ir/Ax 
in  each  direction),  and  k  is  the  magnitude  of  the  vector  wavenumber  k.  In  some  cases,  the  medium 
properties  c(x)  and  p(x)  were  also  smoothed  by  windowing  in  the  spatial-frequency  domain  using 
Eq.  (16)  with  a  wavenumber  cutoff  of  kmax  =  ft /Ax. 

For  comparison,  wavefields  were  also  computed  using  a  second-order  in  time,  fourth-order  in 
space  finite-difference  method,  described  in  Refs.  [21],  [23],  and  [25]-[27].  As  for  the  fc-space 
computations,  the  incident  wave  was  specified  by  a  single  initial  condition  rather  than  updated  at 
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each  time  step.  Periodic  boundary  conditions  were  applied  on  all  sides  of  the  grid.  Time  steps  were 
determined  using  a  CFL  number  of  0.25,  which  is  a  natural  choice  for  this  finite-difference  method 
[26].  As  in  Refs.  [23]  and  [28],  computations  were  performed  at  each  time  step  only  on  portions 
of  the  grid  where  the  wavefields  were  nonzero;  this  reduces  the  required  computation  time  for  the 
finite-difference  method  by  about  one  half. 

To  test  the  A;-space  and  finite-difference  methods  quantitatively,  benchmark  computations  were 
performed  using  an  exact  series  solution  for  the  scattering  of  a  plane  wave  by  a  fluid  cylinder 
[29].  The  sampling  rate  and  waveform  shape  were  chosen  to  match  the  time-domain  simulation 
data  for  the  case  of  interest.  The  pressure  field  was  then  computed  for  each  frequency  component 
with  relative  magnitude  within  60  dB  of  the  peak  magnitude.  Each  single-frequency  computation 
truncated  the  series  at  the  term  having  a  relative  contribution  less  than  10-12  times  the  sum  of  all 
terms.  The  frequency-domain  scattered  fields  were  then  inverted  by  FFT  to  obtain  numerically 
exact  solutions  for  the  time-domain  pressure  fields  at  the  simulated  measurement  points.  An  exact 
time-domain  solution  for  scattering  from  a  fluid  sphere  was  also  obtained  using  an  analogous 
approach. 

Benchmark  studies  of  accuracy  were  performed  using  a  cylinder  with  radius  2.0  mm  and  acous¬ 
tic  properties  of  human  fat,  and  a  background  medium  with  acoustic  properties  of  water  at  body 
temperature.  Rationale  for  use  of  these  values  is  discussed  in  Ref.  [23].  The  cylinder  had  a  sound 
speed  of  1.478  mm//zs  and  a  density  of  0.950  g/cm3,  while  the  background  medium  had  a  sound 
speed  of  1.524  mm//xs  and  a  density  of  0.993  g/cm3.  The  scattering  geometry  was  as  shown  in 
Fig.  1.  The  incident  pulse  was  a  plane  wave  with  Gaussian  temporal  characteristics,  a  temporal 
Gaussian  parameter  a  =  0.25  /zs,  and  a  central  starting  position  of  x  =  -4.5  mm  at  time  zero. 
For  this  pulse,  a  nominal  maximum  frequency  is  4.43  MHz,  corresponding  to  the  spectral  point 
40  dB  down  from  the  center  frequency  (for  the  benchmark  problem,  this  frequency  corresponds 
to  a  minimum  wavelength  of  0.33  mm).  The  A;-space,  leapfrog  pseudospectral,  finite-difference, 
and  exact  methods  described  above  were  used  to  compute  time  histories  of  the  total  pressure  field 
at  128  equally-spaced  “measurement”  points  spanning  a  circle  of  radius  2.5  mm  concentric  to 
the  cylinder.  The  pressure  was  interpolated  using  a  two-dimensional  lowpass  interpolation  filter 
implemented  by  the  formula  [30] 
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Pinterp  (x,  y) 


^  ^  sin(7r(x  -  Xj)J_Ax)  I0[/3(l  -  [(s  -  Xj)/(mAx)]2)1/2] 

\t  \  7r(x  —  Xi)/ Ax  h[P] 

..  sin(7r(y  -  Vi) /Ax)  I0[£(l  -  [(?/  -  yi)l{mAx))2)ll2}  w  , 

X  — ~7~ - _ .  \  /  a  __ - TToi  x  P\xt,yt) 


x  —  mAx  <xt<x  +  mAx, 
y  —  mAx  <yi  <y  +  mAx, 

where  J0  is  the  zeroth-order  modified  Bessel  function  of  the  first  kind  and  /3  is  the  Kaiser  window 
coefficient,  taken  here  to  be  7.0.  This  choice  of  fi  provides  a  filter  with  flat  response  up  to  about 
0.6  kmax  and  sidelobes  at  the  —70  dB  level. 

The  domain  size  for  each  k- space,  pseudospectral,  and  finite-difference  computation  employ¬ 
ing  this  cylinder  was  18  x  18  mm2. 

Further  studies  of  accuracy  were  performed  using  a  cylinder  of  radius  10  mm.  Other  parameters 
were  as  described  above  for  the  small  problem,  except  that  the  radius  of  the  measurement  circle  was 
12.5  mm  and  the  starting  position  of  the  wavefront  was  x  =  —14.5  mm.  The  k- space  method  was 
employed  to  compute  two  cases  corresponding  to  unsmoothed  and  smoothed  contrast  functions, 
using  a  spatial  step  of  four  points  per  minimum  wavelength  and  a  CFL  number  of  0.5.  In  each 
&-space  computation  for  this  cylinder,  the  domain  size  employed  was  72  x  72  mm2.  The  finite- 
difference  method  was  employed  to  compute  a  single  case,  using  a  spatial  step  of  fourteen  points 
per  minimum  wavelength,  a  CFL  number  of  0.25,  and  a  domain  size  of  72  mm  x  60  mm. 

To  evaluate  the  relative  accuracy  and  efficiency  of  the  A;-space  and  finite-difference  methods 
for  a  high-contrast  scatterer,  computations  were  also  performed  using  a  cylinder  of  radius  2.0  mm 
with  the  sound  speed  and  density  of  human  bone.  The  values  employed  were  a  sound  speed 
of  3.54  mm///.s  and  a  density  of  1.99  g/cm3,  as  in  Ref.  [28].  The  incident  pulse,  receiver,  and 
computational  domain  characteristics  were  identical  to  those  for  the  2.0  mm  “fat”  cylinder  case 
described  above. 

In  all  of  the  above  accuracy  tests,  a  quantitative  measure  of  the  accuracy  was  obtained  us¬ 
ing  the  time-domain  L 2  error  of  each  numerically  computed  pressure  field  pn um(x,  t)  versus  the 
corresponding  exact  series  solution  pexact(x,  t).  This  quantity  has  the  definition 

_  |Pnum  (*r,  ~  P exact  (xr  >  0 11  1201 

6_  bexact(xr,*)||  ’  1 
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where  |p(xr,  t)  ||  is  the  L2  norm  [31]  of  a  matrix  composed  of  the  time-domain  signal  p(x,  t)  for  all 
receiver  points  xr  and  all  time  samples  computed.  Eq.  (20)  represents  an  accuracy  criterion  that 
is  much  stricter  than  more  general  criteria,  such  as  comparison  of  the  rms  waveform  amplitude  or 
the  amplitude  and  phase  at  the  center  frequency.  To  achieve  a  low  L2  error  by  the  definition  of 
Eq.  (20),  both  the  waveform  amplitude  and  phase  must  be  accurately  computed  for  all  significant 
frequency  components  of  the  field. 

The  use  of  the  present  fc-space  method  in  a  more  realistic  two-dimensional  simulation  of  ul¬ 
trasonic  propagation  was  also  tested.  For  this  purpose,  a  cross-sectional  tissue  map  of  the  human 
chest  wall  [28]  was  used  as  the  simulated  medium.  A  pulse  center  frequency  of  3.0  MHz  was 
employed  together  with  a  temporal  Gaussian  parameter  of  0.3127  ps\  these  parameters  correspond 
to  the  highest  center  frequency  employed  in  the  simulation  study  reported  in  Ref.  [28].  The  corre¬ 
sponding  nominal  minimum  wavelength  is  0.34  mm.  The  &-space  computation  employed  4  points 
per  minimum  wavelength,  a  CFL  number  of  0.5,  and  a  grid  size  of  54.9  x  54.9  mm2.  The  finite- 
difference  computation  employed  14  points  per  minimum  wavelength,  a  CFL  number  of  0.25,  and 
a  grid  size  of  38.5  x  29.7  mm2.  As  in  Refs.  [23]  and  [28],  periodic  boundary  conditions  were  ap¬ 
plied  on  the  sides  perpendicular  to  the  wavefront,  while  first-order  radiation  boundary  conditions 
[23]  were  applied  on  the  sides  parallel  to  the  wavefront. 

Finally,  to  illustrate  the  efficiency  and  accuracy  of  the  present  k- space  method  for  three- 
dimensional  computations,  scattering  from  a  penetrable  sphere  with  acoustic  properties  of  hu¬ 
man  muscle  (speed  1.547  mm//xs,  density  1.090  g/cm3  [23])  was  computed.  The  sphere  radius 
was  1.5  mm;  time-domain  pressure  waveforms  were  recorded  at  128  equally-spaced  measurement 
points  on  the  sphere  surface  (in  the  (j)  =  0  plane).  The  computation  employed  an  incident  pulse 
identical  to  that  for  the  cylinder  simulations  described  above,  a  spatial  step  of  four  points  per  min¬ 
imum  wavelength,  and  a  CFL  number  of  0.5.  The  total  pressure  wavefield  was  computed  for  a 
time  duration  of  7.3  ps  on  a  three-dimensional  grid  of  dimensions  10.66  x  10.66  x  10.66  mm3. 
The  accuracy  of  this  computation  was  assessed  by  evaluating  the  L2  error  between  the  fc-space  and 
exact  solutions  using  Eq.  (20). 
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IV.  Numerical  Results 


An  example  &-space  computation,  performed  using  the  2.0  mm  cylinder  with  acoustic  properties  of 
human  fat,  is  illustrated  in  Fig.  1.  The  cylinder  is  also  sketched  in  each  panel.  For  the  computation 
shown,  smoothed  sound  speed  and  density  functions  were  obtained  by  filtering  the  analytic  spatial 
Fourier  transform  of  the  cylinder  using  Eq.  (16).  The  time  history  of  the  total  wavefield  is  shown 
as  computed  by  the  A; -space  method  for  a  spatial  step  size  of  four  points  per  minimum  wavelength 
and  a  CFL  number  of  0.5.  Details  visible  include  a  scattered  wave  from  the  edge  nearest  the  initial 
wavefront  (c),  weak  focusing  near  the  trailing  edge  of  the  cylinder  (e),  scattering  from  the  trailing 
edge  [(f)— (i)] ,  and  low-level  multiple  scattering  [(g)— (h)]. 

Results  of  accuracy  benchmarks  for  the  k- space  and  leapfrog  pseudospectral  methods  described 
above  are  shown  in  Fig.  2.  Each  of  these  computations  was  made  using  the  2.0  mm  radius  cylinder 
described  above  and  a  spatial  step  size  of  four  points  per  maximum  wavelength.  The  results  show 
that  the  ifc-space  method  employing  the  k-t  space  propagator  of  Eq.  (9)  provides  much  higher 
accuracy  than  the  pseudospectral  method  employing  the  leapfrog  propagator  of  Eq.  (8).  The  two 
methods  provide  equivalent  results  for  very  small  time  steps  (CFL  numbers  less  than  about  0. 1),  but 
the  A-space  method  maintains  its  highest  accuracy  up  to  a  CFL  number  of  about  0.4.  In  contrast, 
the  pseudospectral  method  rapidly  increases  in  error  for  CFL  numbers  above  0.1. 

Error  results  for  the  pseudospectral  computations  shown  in  Fig.  2  are  not  given  for  CFL  num¬ 
bers  above  0.6  because  the  computation  was  unstable  for  higher  CFL  numbers.  (That  is,  computed 
fields  incurred  spurious  exponential  growth,  resulting  in  numerical  overflow.)  This  observation  of 
instability  is  consistent  with  the  linear  stability  limit  of  0.6366  given  by  Eq.  (12)  for  this  case.  The 
A;-space  method  did  not  incur  any  numerical  instability  for  the  range  of  CFL  numbers  investigated, 
so  that  the  method  is  seen  to  be  unconditionally  stable  as  predicted  for  c(x)  <  cq.  However,  the 
error  of  this  method  grows  as  the  CFL  number  approaches  and  exceeds  unity,  consistent  with  the 
Nyquist  sampling  criterion  given  by  Eq.  (14). 

Pseudospectral  methods  employing  higher-order  time  integration  achieve  higher  accuracy  than 
the  leapfrog  iteration  used  as  a  comparison  here.  However,  tests  of  the  present  fc-space  method  and 
a  pseudospectral  method  employing  fourth-order  Adams-Bashforth  time  integration  have  shown 
trends  similar  to  that  seen  in  Fig.  2  [32].  Specifically,  for  weakly-scattering  media,  the  /c- space 
method  yields  similar  accuracy  for  time  steps  two  to  three  times  larger  than  those  required  by  the 
higher-order  pseudospectral  method  described  in  Ref.  [13]. 
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The  relative  accuracy  of  the  fc-space  method  and  the  2-4  finite  difference  method  are  compared 
in  Fig.  3  as  a  function  of  the  spatial  step  size.  For  these  computations,  the  CFL  number  of  the 
&-space  computations  was  held  constant  at  0.5,  consistent  with  the  CFL-accuracy  relationship 
shown  in  Fig.  2,  while  the  CFL  number  of  the  finite-difference  computations  was  held  at  0.25 
[26].  Both  methods  achieve  high  accuracy  for  finer  grid  spacings;  however,  the  fc-space  method 
achieves  higher  accuracy  for  much  larger  spatial  step  sizes.  The  L2  error  drops  below  0.05  for 
fc-space  computations  employing  only  three  points  per  minimum  wavelength,  while  achievement 
of  the  same  accuracy  criterion  requires  14  points  per  minimum  wavelength  for  the  finite-difference 
computations.  This  difference  suggests  that  storage  requirements  for  k- space  computations  can  be 
much  smaller  than  those  for  finite-difference  computations  of  comparable  accuracy:  on  the  order 
of  12  times  smaller  for  two-dimensional  computations  and  43  times  smaller  for  three-dimensional 
computations. 

Visual  comparison  of  simulated  waveforms  for  the  2.0  mm  radius  cylinder  is  shown  in  Fig.  4. 
Waveforms  in  this  figure  are  those  computed  using  the  k- space  (four  points  per  minimum  wave¬ 
length,  CFL  number  0.5,  with  both  unsmoothed  and  smoothed  contrast  functions),  finite-difference 
time-domain  (14  points  per  minimum  wavelength,  CFL  number  0.25),  and  exact  methods.  The  k- 
space  solution  for  the  unsmoothed  cylinder  shows  a  small  time-domain  L2  error  (0.0243),  but  also 
exhibits  spurious  waves  (nearly  60  dB  down  from  the  peak  pressure  amplitude)  between  the  two 
main  arrivals.  These  spurious  waves  are  removed  by  use  of  the  /c- space  method  with  smoothed 
medium  parameters  [i.e.,  p(x)  and  c(x)  smoothed  using  Eq.  (16)  with  kmax  =  n/Ax];  the  L2  error 
is  decreased  to  0.0214  by  this  smoothing.  The  finite-difference  result  bears  a  strong  qualitative 
resemblance  to  the  exact  solution,  but  the  larger  L2  error  (0.0454)  indicates  that  phase  errors  have 
been  introduced  by  the  dispersion  inherent  to  the  finite-difference  method.  Computation  times  [33] 
were  2.31  minutes  for  the  k- space  method  and  1.55  hours  for  the  finite-difference  method,  so  that 
the  A: -space  method  yields  greater  accuracy  at  much  less  computational  cost. 

Waveforms  for  the  10  mm  radius  cylinder  are  shown  in  Fig.  5  in  a  format  analogous  to  that  of 
Fig.  4.  These  results  indicate  that,  as  for  the  smaller  cylinder,  smoothing  of  the  contrast  functions 
produces  a  reduction  in  spurious  low-amplitude  waves.  For  this  problem,  unlike  the  2.0  mm  radius 
cylinder  discussed  above,  this  smoothing  slightly  decreases  the  overall  accuracy.  (The  time-domain 
L2  error  is  0.1292  for  the  smoothed  case  versus  0.1288  for  the  unsmoothed  case.)  The  finite- 
difference  solution,  using  14  points  per  wavelength  and  a  CFL  number  of  0.25,  requires  much 
greater  storage  and  computational  time,  and  produces  waveforms  with  poorer  accuracy  (an  L2 
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error  of  0.1794)  than  the  /c-space  method. 

Results  for  the  2  mm  “bone”  cylinder  are  shown  in  Fig.  6.  In  this  case,  the  A: -space  method 
using  a  CFL  number  of  0.5  exhibited  numerical  instability.  This  instability  is  expected,  since 
this  CFL  number  exceeds  the  limit  of  0.2833  set  by  Eq.  (13).  To  obtain  an  appropriate  temporal 
sampling  rate,  the  time  step  was  reduced  in  proportion  to  the  increase  in  cmax,  resulting  in  a  CFL 
number  of  0.2153.  Required  computation  time  for  the  A:-space  method  was  5.34  minutes  [33];  the 
time-domain  L 2  error  was  0.3061  for  the  unsmoothed  case  and  0.2687  for  the  smoothed  case. 

The  finite-difference  method,  employing  14  points  per  wavelength  and  a  CFL  number  of  0.1076 
(also  changed  in  proportion  to  Cm  ax),  achieved  an  L2  error  of  0.0350  in  a  computation  time  of 
3.99  hours  [33].  This  result  indicates  that  finite-difference  methods  can  be  much  more  accurate 
than  A: -space  methods  for  scattering  problems  involving  very  high-contrast  inhomogeneities  such 
as  bone  within  soft  tissue.  However,  the  k- space  solution,  as  seen  in  Fig.  6,  still  shows  good 
qualitative  agreement  with  the  exact  solution. 

The  relative  inaccuracy  of  the  A;-space  method  for  high-contrast  scatterers  may  be  associated 
with  aliasing  effects,  as  suggested  in  Ref.  [5].  That  is,  large  jumps  in  spatial  contrast  functions  are 
associated  with  significant  high-frequency  components  of  the  corresponding  A;-space  spectra.  If  the 
spatial-frequency  range  employed  in  the  /c-space  algorithm  is  not  sufficiently  large,  aliasing  errors 
result.  Low-pass  filtering  of  the  contrast  functions  would  remove  this  aliasing,  but  also  introduces 
additional  errors  because  high  spatial-frequency  components  of  the  scattering  medium  are  lost. 
The  half-band  filtering  employed  here  is  a  compromise  that  greatly  reduces  aliasing  errors  while 
maintaining  some  contributions  from  high  spatial  frequencies  (up  to  the  spatial  Nyquist  rate). 

Computational  results  for  a  large-scale  two-dimensional  tissue  model  are  shown  in  Fig.  7. 
Waveforms  computed  by  the  A;- space  (four  points  per  minimum  wavelength,  CFL  number  0.5,  no 
smoothing)  and  the  finite-difference  (ten  points  per  minimum  wavelength,  CFL  number  0.25)  were 
recorded  at  130-element  apertures  composed  of  simulated  point  receivers  separated  by  a  pitch  of 
0.21  mm.  The  results  produced  by  the  finite-difference  method  and  the  A;-space  method  are  visually 
indistinguishable.  However,  despite  the  reduced  grid  size  and  limited  computations  employed  for 
the  finite-difference  method,  the  A;-space  method  was  more  efficient  by  about  a  factor  of  four;  the 
required  CPU  time  for  the  A;-space  method  was  0.90  CPU  hours,  while  the  corresponding  time  for 
the  finite-difference  time-domain  method  was  4.58  CPU  hours  [33].  This  discrepancy  in  efficiency 
is  even  more  impressive  when  note  is  made  that  the  A;-space  method  using  4  points  per  minimum 
wavelength  provides  significantly  higher  accuracy  than  the  finite-difference  method  using  14  points 
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per  minimum  wavelength  (as  illustrated  in  Fig.  3).  Thus,  the  present  fc-space  method  is  suggested 
to  be  an  appropriate  replacement  for  finite-difference  methods  previously  employed  to  compute 
propagation  through  large-scale  soft-tissue  models  [23]-[28]. 

Results  of  the  example  three-dimensional  computation  are  shown  in  Fig.  8.  Three-dimensional 
isosurface  renderings  of  the  total  pressure  wavefield  are  shown  at  three  instants  separated  by 
0.79  ^s.  For  the  three-dimensional  computation,  the  total  computation  time  required  was 
1.51  hours  [33].  The  L2  error  of  the  computed  waveforms,  relative  to  the  exact  time-domain 
solution  for  scattering  from  a  sphere  [29],  was  0.0186. 
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V.  Extensions  to  the  £>Space  Method 


The  present  method  can  be  extended  in  a  number  of  ways  to  increase  its  range  of  applicability  in 
computations  of  ultrasound-tissue  interactions. 

Absorption  effects  could  be  added  to  the  present  algorithm  in  several  ways.  The  most  straight¬ 
forward  method  for  including  absorption  is  to  include  an  ad  hoc  damping  term  proportional  to 
dfa/dt  in  Eq.  (2)  [3]-[5].  This  approach  yields  absorption  coefficients  roughly  independent  of 
the  frequency.  Similarly,  inclusion  of  a  damping  term  proportional  to  dzfs/dtz  (a  thermovis- 
cous  approximation)  would  lead  to  absorption  roughly  proportional  to  the  frequency  squared  [34]. 
However,  neither  of  these  approaches  has  a  rigorous  justification  for  use  in  models  of  ultrasound 
propagation  in  biological  tissue. 

A  physically  justifiable  approach  for  inclusion  of  absorption  in  the  present  algorithm  is  to 
consider  absorption  associated  with  multiple  relaxation  processes.  The  theoretical  basis  for  this 
approach  is  presented  in  Ref.  [35];  one  implementation  of  this  method  in  a  finite-difference  time- 
domain  algorithm  is  given  in  Ref.  [36].  Since  multiple  relaxation  processes  can  lead  to  a  variety  of 
frequency-dependent  absorption  characteristics,  this  approach  provides  a  possibility  of  modeling 
realistic  frequency-dependent  attenuation  in  tissue  without  introduction  of  nonphysical  dispersion 
or  violation  of  causality.  Following  the  methods  presented  in  Ref.  [36],  absorption  due  to  multiple 
relaxation  processes  can  be  implemented  in  a  computationally  efficient  form.  Possible  alternatives 
include  the  time-causal  power  law  absorption  formulation  of  Ref.  [37]. 

Another  possible  extension  to  the  present  method  is  to  incorporate  the  full  elastic  wave  propa¬ 
gation  equations.  This  extension  would  account  for  shear  wave  propagation,  which  may  substan¬ 
tially  affect  results  for  propagation  models  including  bone  and  other  calcified  tissue.  By  applying 
methods  similar  to  those  outlined  in  Ref.  [7]  to  the  algorithm  described  above,  a  full  elastic  k- space 
method  incorporating  Fourier-space  evaluation  of  spatial  derivatives  and  a  k-t  space  propagator 
could  be  derived.  Such  a  method  would,  as  in  Ref.  [7],  include  separate  k-t  space  propagators  for 
compressional  and  shear  waves. 

Boundary  conditions  of  k- space  and  pseudospectral  methods  are  inherently  periodic,  so  that 
simple  radiation  boundary  conditions  cannot  be  straightforwardly  implemented.  One  option  for  ab¬ 
sorbing  boundary  conditions  is  to  include  tapered  (artificial)  absorption  functions  at  each  boundary 
[38].  The  technique  of  perfectly  matched  layers  (PML)  [39]  can  provide  true  radiation  boundary 
conditions;  however,  present  PML  implementations  are  not  applicable  to  the  second-order  wave 
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equation  employed  here.  Combination  of  a  k- space  method  with  PML  boundary  conditions  may 
require  derivation  of  a  new  k-t  space  time  integrator  for  the  first-order  wave  propagation  equations. 

The  present  derivation  was  based  on  the  linear  (small-amplitude)  acoustic  propagation  equa¬ 
tions.  The  6-space  method  could  be  easily  extended  to  incorporate  finite-amplitude  acoustic  ef¬ 
fects.  For  example,  the  nonlinear  terms  of  the  Westervelt  propagation  equation  (used  in  Ref.  [34] 
for  modeling  of  ultrasonic  propagation  in  tissue)  could  be  included  as  effective  source  terms  addi¬ 
tional  to  the  effective  sources  v  and  q  defined  above.  The  numerical  results  obtained  above  suggest 
that  the  6-space  method  is  most  accurate  when  the  effective  source  terms  are  fairly  small;  thus,  a 
nonlinear  extension  to  the  6 -space  method  should  be  highly  accurate  for  weakly  nonlinear  effects. 

Computation  times  for  the  k- space  method  can  easily  be  reduced  by  parallelization.  The  pri¬ 
mary  computational  burden  of  the  method  is  incurred  in  the  multidimensional  fast  Fourier  trans¬ 
forms  (FFT)  taken  at  each  time  step.  Since  FFT’s  can  be  efficiently  executed  on  parallel  processors 
[24, 40],  the  present  6-space  method  should  scale  efficiently  to  large  problems  that  require  parallel 
processing. 
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VI.  Conclusions 


A  simplified  derivation  of  the  A;-space  method  for  computation  of  ultrasonic  wave  propagation 
has  been  presented.  The  method  efficiently  accounts  for  sound  speed  and  density  variations,  and 
can  be  extended  to  include  realistic  absorption  effects  and  absorbing  boundary  conditions.  Three- 
dimensional  computations  can  also  be  performed  without  change  to  the  algorithm  as  derived  here. 

Analytic  and  numerical  results  have  shown  that  the  present  A;-space  method  provides  superior 
stability  and  accuracy  over  both  a  similar  leapfrog  pseudospectral  method  and  a  fourth-order  space, 
second-order  time,  finite-difference  method.  This  improved  accuracy  allows  larger  spatial  and 
time  steps  to  be  employed,  so  that  large-scale  multidimensional  computations  are  more  feasible. 
Computations  using  a  realistic  two-dimensional  tissue  model  support  the  conclusion  that  the  k- 
space  method  provides  high  accuracy  and  low  computational  cost  for  large-scale  computations. 

The  results  also  indicate  that  care  should  be  taken  when  choosing  and  implementing  a  forward 
solver  for  a  particular  scattering  problem.  For  instance,  in  the  present  k- space  method,  one  can 
suppress  spurious  waves  by  smoothing  sound  speed  and  density  variations;  however,  this  smooth¬ 
ing  does  not  decrease  the  time-domain  L2  error  in  some  cases.  Likewise,  the  finite-difference 
time-domain  method  employed  here  is  less  accurate  than  the  A;-space  method  in  most  cases  ex¬ 
amined  here,  but  achieved  higher  accuracy  for  a  test  case  with  a  bone-like  scatterer.  In  general, 
the  fc-space  method  proposed  here  should  be  most  applicable  to  large-scale  scattering  problems 
involving  low-contrast  inhomogeneities  such  as  soft  tissue  structures. 
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Figure  1:  Time  history  of  total  acoustic  pressure  computed  by  the  k- space  method  for  a  cylinder 
of  2.0  mm  radius  and  fat-mimicking  acoustic  properties.  The  cylinder  is  sketched  as  a  light  gray 
region.  The  first  panel  shows  the  wavefield  impinging  on  the  cylinder  at  time  t  =  0.98  \i s  and 
subsequent  panels  (progressing  from  left  to  right  and  top  to  bottom)  show  the  total  wavefield  at 
intervals  of  0.98  //s.  The  acoustic  pressure  is  plotted  in  all  panels  using  a  bipolar  logarithmic  scale 
with  a  60  dB  dynamic  range. 


Figure  2:  Time-domain  comparison  of  accuracy  for  the  fc-space  and  leapfrog  pseudospectral  meth¬ 
ods  as  a  function  of  CFL  number.  Each  test  used  the  “fat”  cylinder  of  2.0  mm  radius  and  a  spatial 
step  size  of  four  points  per  minimum  wavelength. 
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Figure  3:  Time-domain  comparison  of  accuracy  for  the  k- space  and  2-4  finite-difference  time- 
domain  methods  as  a  function  of  the  spatial  step  size  in  points  per  minimum  wavelength  (PPW). 
Each  test  used  the  “fat”  cylinder  of  2.0  mm  radius.  CFL  numbers  were  0.5  for  the  fc-space  method 
and  0.25  for  the  finite-difference  time-domain  method. 
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Figure  4:  Computed  waveforms  for  the  “fat”  cylinder  at  a  radius  of  2.5  mm  for  a  cylinder  of  radius 
2.0  mm  and  a  pulse  center  frequency  of  2.5  MHz.  The  acoustic  pressure  is  shown  on  a  bipolar 
logarithmic  scale  with  60  dB  dynamic  range.  The  horizontal  range  of  each  plot  is  360  degrees, 
covering  the  entire  measurement  circle  starting  with  angle  0  (forward  propagation).  The  vertical 
range  of  each  panel  corresponds  to  a  temporal  duration  of  9  (jls,  with  t  =  0  at  the  top  of  each 
plot,  (a)  Unsmoothed  object;  k- space  solution  with  four  points  per  minimum  wavelength,  L2  error 
0.0243.  (b)  Smoothed  object;  fc-space  solution  with  four  points  per  minimum  wavelength,  L2  error 
0.0214.  (c)  Finite-difference  solution  with  14  points  per  minimum  wavelength,  L2  error  0.0454. 
(d)  Exact  solution. 


Figure  5:  Computed  waveforms  at  a  radius  of  12.5  mm  for  a  “fat”  cylinder  of  radius  10.0  mm 
and  a  pulse  center  frequency  of  2.5  MHz.  The  acoustic  pressure  is  shown  in  each  panel  using  a 
bipolar  logarithmic  scale  with  a  60  dB  dynamic  range.  The  horizontal  range  of  each  panel  is  360 
degrees  and  the  vertical  range  is  33  //s.  (a)  Unsmoothed  object;  fc-space  solution,  L2  error  0.1288. 
(b)  Smoothed  object;  fc-space  solution,  L2  error  0.1292.  (c)  Finite-difference  solution,  L2  error 
0.1794.  (d)  Exact  solution. 
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Figure  6:  Computed  pressure  waveforms  at  a  receiver  radius  of  2.5  mm  for  a  “bone”  cylinder  of 
radius  2.0  mm  and  a  pulse  center  frequency  of  2.5  MHz.  The  format  is  the  same  as  in  Fig.  4.  (a) 
Unsmoothed  object;  &-space  solution,  L2  error  0.3061.  (b)  Smoothed  object;  k- space  solution,  L2 
error  0.2687.  (c)  Finite-difference  solution,  L2  error  0.0380.  (d)  Exact  solution. 


Figure  7:  Comparison  of  /c-space  and  finite-difference  methods  for  a  tissue  cross-sectional  model, 
(a)  Chest  wall  cross  section  (taken  from  Ref.  [24]),  with  black  indicating  connective  tissue,  dark 
gray  indicating  muscle,  and  light  gray  indicating  fat.  The  region  is  33.5  mm  wide  and  17.2  mm 
high,  (b)  Transmitted  waveforms  computed  by  the  k- space  method  using  four  points  per  minimum 
wavelength  and  a  CFL  number  of  0.5,  shown  on  a  bipolar  linear  gray  scale  with  white  indicating 
maximum  positive  pressure  and  black  indicating  maximum  negative  pressure.  The  horizontal  range 
shown  is  27.3  mm  and  is  shown  to  the  same  scale  as  in  (a).  The  vertical  range  is  3.29  fj,s.  (c) 
Transmitted  waveforms  computed  by  the  finite-difference  time-domain  method  using  10  points 
per  minimum  wavelength  and  a  CFL  number  of  0.25,  shown  using  the  same  format  as  in  (b). 
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Figure  8:  Isosurface  renderings  of  the  total  (logarithmically  scaled)  pressure  wavefield  associated 
with  scattering  from  a  “muscle”  sphere  of  radius  1.5  mm.  Incident  pulse  parameters  were  the 
same  as  in  Figs.  4-6.  Panels  (a)-(d)  show  the  wavefield  at  four  instants  separated  by  0.79  //s.  The 
view  shown  is  such  that  the  incident  wave  is  traveling  into  the  page,  so  that  the  visible  wavefield 
includes  the  backscattered  component.  The  lowest-amplitude  isosurface  shown  is  67.5  dB  down 
from  the  incident-wave  amplitude.  Each  panel  shows  a  rendering  of  the  entire  computational 
domain  (10.66  mm  on  each  side).  In  panel  (a);  the  incident  wavefront  is  just  impinging  on  the 
sphere;  in  panel  (d),  the  scattered  wavefront  has  just  passed  the  computational  boundary. 
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2aBB6.  A  new  it -space  method  for  simulation  of  ultrasonic 
propagation  in  tissue.  T.  Douglas  Mast  (Appl.  Res.  Lab.,  Penn  State 
Univ.,  University  Park,  PA  16802,  mast@sabine.acs.psu.edu),  D.-L. 
Donald  Liu  (Siemens  Medical  Systems,  Issaquah,  WA  98027),  Laurent 
P.  Souriau,  Adrian  I.  Nachman,  and  Robert  C.  Waag  (Univ.  of  Rochester, 
Rochester,  NY  14642) 

A  new  k-space  method  for  large-scale  computations  of  ultrasonic 
propagation  is  presented.  In  the  new  method,  spatial  derivatives  from  the 
second-order  acoustic  wave  equation  for  inhomogeneous  media  are  evalu¬ 
ated  by  Fourier  transformation.  Solutions  are  advanced  in  time  using  a 
k-t  space  Green’s  function.  Computational  results  indicate  that  the  new 
method  shares  advantages  of  both  past  fc-space  and  pseudospectral  meth¬ 
ods.  For  scatterers  with  properties  similar  to  soft  tissue,  the  fc- space 
method  provides  much  higher  accuracy  and  lower  computational  cost  than 
a  2-4  finite-difference  time  domain  method.  The  k-space  method  also 
allows  high  accuracy  to  be  obtained  for  time  steps  much  larger  than  those 
required  by  a  leapfrog  pseudospectral  method.  The  low  dispersion  inherent 
to  the  & -space  method  is  illustrated  by  large-scale  quasi-one-dimensional 
computations,  in  which  pulse  waveforms  incur  negligible  shape  change  for 
propagation  distances  as  large  as  1000  wavelengths.  Example  applications 
of  the  k-space  method  are  demonstrated,  including  simulation  of  propaga¬ 
tion  through  a  large-scale  tissue  cross-sectional  model  and  incorporation 
of  a  t-space  solver  into  a  nonlinear  inverse  scattering  method  employing 
eigenfunctions  of  the  far-field  scattering  operator. 


TIME-DOMAIN  INVERSE  SCATTERING  FOR 
QUANTITATIVE  ULTRASONIC  MAMMOGRAPHY 

T.  Douglas  Mast 


Applied  Research  Laboratory 
The  Pennsylvania  State  University 
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A  new  method  for  ultrasonic  mammography  is  presented.  This  method  provides 
quantitative  tomographic  images  of  inhomogeneous  tissue  using  time-domain  scattering 
measurements  made  on  a  surrounding  surface  (for  example,  on  a  circle  for  images  of  a 
two-dimensional  breast  cross  section).  High-resolution,  quantitative  images  of  tissue  are 
reconstructed  using  coherent  combination  of  far-field  scattered  ultrasound  waveforms, 
delayed  and  summed  to  focus  at  each  image  point.  The  focused  image  is  a  reconstruction 
of  the  spatially-dependent  sound  speed  variation,  and  is  equivalent  to  a  wideband  filtered 
backpropagation  reconstruction  weighted  by  the  spectrum  of  the  incident  wave.  The 
resulting  images  are  higher  in  quality  than  frequency-domain  quantitative  reconstructions 
and  contain  more  diagnostic  information  than  conventional  B-scans. 

Rigorous  testing  of  the  new  imaging  method  is  carried  out  using  simulated  ultrasonic 
propagation  through  breast  tissue.  Breast  tissue  models  are  obtained  both  from 
segmentation  of  stained  cross  sections  and  from  analysis  of  high-resolution  three- 
dimensional  data  from  the  Visible  Woman  project.  Computations  of  ultrasonic 
propagation  are  performed  using  a  new  k- space  method,  in  which  the  spatial  differential 
equations  are  solved  by  Fourier  transformation  and  temporal  iteration  is  performed  using 
a  k-t  space  propagator.  Numerical  results  indicate  that  this  method  is  highly  accurate  for 
large-scale  soft-tissue  computations,  with  much  greater  efficiency  than  that  of  competing 
methods.  Thus,  the  k-space  method  is  particularly  appropriate  for  large-scale  two- 
dimensional  and  three-dimensional  computations  of  propagation  through  breast  tissue. 

Quantitative  images,  obtained  using  synthetic  data  for  two-dimensional  and 
three-dimensional  scattering  of  wideband  pulses  as  well  as  measured  scattering  data  from 
a  2048-element  ring  transducer,  confirm  that  the  time-domain  reconstruction  method 
provides  superior  image  quality  for  objects  of  size  and  contrast  relevant  to  ultrasonic 
mammography.  The  new  method  can  also  be  extended  to  incorporate  available 
image-enhancement  techniques,  such  as  time-gain  compensation  to  correct  for  medium 
absorption  and  aberration  correction  methods  to  reduce  error  associated  with  weak 
scattering  approximations. 
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New  methods  for  breast  cancer  detection  and  characterization  are  the  focus  of  this  project, 
supported  by  the  US  Army  Medical  Research  and  Materiel  Command.  Major  goals  are 
to  establish  a  new  high-resolution,  quantitative  ultrasonic  imaging  method  and  to  test  this 
method  using  simulated  propagation  of  ultrasonic  pulses  through  accurately  detailed  breast 
tissue  models. 

Detailed  breast  tissue  models  have  been  obtained  both  from  stained  breast  cross  sections 
and  from  analysis  of  high-resolution  three-dimensional  data  from  the  Visible  Woman 
project.  To  accurately  compute  ultrasonic  propagation  through  these  tissue  models,  a  new 
it- space  method  for  ultrasound  simulation  has  been  developed.  The  &- space  method  is 
more  accurate,  more  efficient,  and  requires  less  storage  than  alternative  methods,  and  is 
thus  ideal  for  computation  of  large-scale  2D  and  3D  ultrasonic  propagation  in  breast  tissue. 
A  new  time-domain  ultrasonic  mammography  method  provides  quantitative  images  of 
inhomogeneous  media  including  breast  tissue.  High-resolution  maps  of  the  tissue  sound 
speed  are  obtained  from  processing  of  measured  ultrasonic  scattering.  Unlike  previous 
frequency-domain  inverse  scattering  methods,  the  entire  signal  bandwidth  is  used,  so  that 
reconstructed  images  have  higher  point  resolution  (ability  to  detect  small  structures  such 
as  microcalcifications)  and  contrast  resolution  (ability  to  distinguish  subtle  differences 
between  tissue  structures).  The  new  method  employs  a  straightforward  time-domain 
reconstruction  algorithm,  similar  to  synthetic-aperture  methods  used  by  current  clinical 
scanners,  but  provides  much  more  diagnostic  information  than  current  B-scan  devices. 
The  high  efficiency  of  the  reconstruction  algorithm  makes  the  new  method  particularly 
well-suited  for  three-dimensional  quantitative  ultrasonic  mammography. 

Quantitative  images,  obtained  both  from  synthetic  and  measured  ultrasound  data,  confirm 
that  the  new  imaging  method  provides  superior  image  quality  and  accurate  quantitative 
information.  After  further  development  and  clinical  implementation,  the  new  ultrasonic 
mammography  method  is  expected  to  become  competitive  with  magnetic  resonance 
imaging  and  x-ray  computed  tomography  as  a  tool  for  breast  cancer  detection  and 
characterization,  while  maintaining  inherent  advantages  of  ultrasound  such  as  lower  cost, 
ability  to  characterize  cystic  and  solid  lesions,  and  safe,  nonionizing  radiation. 
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